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ABSTRACT 


Three  different  approaches  to  the  problem  of  imple- 
menting a  reduced-order,  sub-optimal  Kalman  filter  for  a 
discrete,  linear  stochastic  process,  with  time-invariant 
dynamics,  are  presented. 

A  first  method,  A,  is  based  upon  the  partitioning  of 
the  system  dynamics.   A  second  method,  B,  is  implemented 
using  matrix  pseudo-inversion  and  a  third  method,  C,  is 
based  upon  reduction  of  the  original  process  to  one  of 
lower  order  using  the  dominant  roots  of  the  system.   An 
expression  for  the  performance  degradation  in  method  A  is 
derived.   In  method  B,  expressions  for  the  sub-optimal 
estimation  error,  and  sub-optimal  variance  of  estimation 
error  are  derived. 

The  several  methods  are  applied  to  a  fourth-order 
process  for  illustration. 
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I.  INTRODUCTION 

In  any  modern  control  problem  involving  filtering  of 
discrete  or  continuous  stochastic  linear  processes,  it  is 
clear  that  the  computational  requirements  associated  with 
the  implementation  of  a  Kalman  filter  increase  consider- 
ably as  the  dimension  of  the  system  increases.   A  simpler 
(reduced)  configuration  of  the  filter  with  sub-optimal 
performance  may  be  acceptable,  mainly  if  some  preselected 
r  states  (where  r-=n ,  and  n  represents  the  order  of  the 
original  system)  are  of  primary  interest  in  the  estimation 
process  and  if  the  performance  degradation  incurred  when 
using  the  reduced  model  is  not  severe. 

With  this  idea  of  implementing  a  reduced,  discrete, 
sub-optimal  filter,  based  upon  the  optimal  Kalman  estima- 
tion theory  (Ref. 1.),  a  library  search  was  done,  yielding 
two  papers  of  particular  interest  which  present  different 
approaches  to  the  problem  of  reduced  estimation  for  con- 
tinuous dynamic  processes.   Meditch  (Ref. 3)  in  1964,  de- 
veloped a  theory  of  sub-optimal  filtering  for  continuous 
dynamic  processes,  in  which  the  process  vector  is  divided 
by  means  of  matrix  pseudo  inversion,  with  separate  filters 
designed  for  each  part.   The  sub-optimal  estimates  are  ob- 
tained by  combining  the  several  filter  outputs.   As 
Meditch  points  out  in  his  paper,  an  important  feature  of 
his  continuous  formulation  is  that  it  provides  a  powerful 
tool  for  preliminary  analysis  of  sub-optimal  fiiters. 


Huddle  and  Wismer  (Ref.4)  in  1968  based  a  sub-optimal 
filter  design  on  a  "primary  system",  which  is  implemented 
by  dividing  the  original  process  into  two  separate  systems, 
using  partitioning  of  the  continuous  dynamics  of  the  sys- 
tem, and  neglecting  the  remaining  "secondary  system". 

Another  successful  technique  for  simplifying  filter 
computations  has  been  investigated  by  Aoki  and  Huddle 
(Ref.5),  wherein  a  procedure  is  developed  to  estimate  the 
state  vector  of  a  discrete  stochastic  system,  when  con- 
straints are  imposed  on  the  number  of  memory  elements  of 
the  estimator. 

This  thesis  presents  both  an  extension  of  two  of  the 
methods  referred  to  above  (Ref.3  and  Ref.4)  to  the  discrete 
case,  and  a  third  approach  to  the  same  problem  based  upon 
the  use  of  the  dominant  roots  of  the  full  process. 

The  study  of  the  three  different  approaches  to  the 
reduced  sub-optimal  filter  problem  is  presented  in  a  co- 
herent form  as  follows.   Chapter  II  describes  the  general 
discrete  linear  stochastic  process.   In  Chapter  III,  the 
work  of  Huddle  and  Wisroer  (Ref.4)  is  extended  to  the  dis- 
crete case,  and  the  performance  degradation  of  the  re- 
duced filter  is  obtained.   An  example  using  a  fourth- 
order  process  is  solved  estimating  only  two  states  of  the 
full  system,  in  order  to  illustrate  the  method.   In  Chap- 
ter IV,  the  Meditch  (Ref.3)  reduced  filter  is  adapted  to 
the  discrete  case,  and  recursive  equations  for  the  sub- 
optimal  estimation  error  and  sub-optimal  covariance  of 


estimation  error  are  derived  and  used  as  a  means  of  com- 
paring the  sub-optimal  with  optimal  performances.   The 
same  fourth-order  system  used  in  Chapter  III,  is  used  here 
as  an  example,  and  again,  only  two  states  are  estimated. 
In  Chapter  V,  another  reduced  filter  is  implemented,  using 
the  dominant  roots  of  the  complete  system  in  order  to  pro- 
duce the  dynamics  of  the  reduced  filter.   Once  more,  the 
technique  is  illustrated  using  the  original  plant  of  the 
previous  examples. 

In  the  three  examples  presented  as  illustrations,  the 
same  complete  process  is  used  in  order  to  compare  the  re- 
sults of  the  three  different  methods.   In  each  one  of  the 
examples,  a  Monte  Carlo  simulation  of  the  full  process  was 
performed,  and  optimal  Kalman  estimation  of  the  whole 
state  vector  was  accomplished  in  order  to  establish  the 
performance  degradation  of  the  reduced  filters. 

Appendix  A  presents  the  basic  principles  of  matrix 
pseudo-inversion;   more  complete  treatments  of  this  sub- 
ject appear  in  Penrose  (Ref.8)  and  Werther  (Ref.9). 

Appendix  B  and  C  include  the  computer  programs  used 
in  solving  the  examples.   These  programs  were  not  intended 
to  be  computationally  efficient  either  in  time,  or  in 
storage. 


II.  THE  DISCRETE  (DISCRETIZED)  LINEAR 
STOCHASTIC  PROCESS  WITH 
TIME- INVARIANT  DYNAMICS 

A  linear,  continuous  stochastic  dynamic  process  may 
be  described  in  general  (Ref.7)  by  the  vector  matrix  dif- 
ferential equation 

x(t)  =  A(t)x(t)  +  B(t)[u(t)+w(t)]  (2.1) 

with  output 

y_(t)  =  H(t)x(t)  (2.2) 

where:  x(t)  is  a  (n  x  1)  vector,  representing  the  state 
variables  of  the  system 

u(t)  is  a  (m  x  1)  vector  (m  n),  representing  the 
deterministic  input  to  the  system 

w(t)  is  a  (m  x  1)  vector  representing  the  stochas- 
tic input  to  the  system 

A(t)  is  a  (n  x  n)  system  matrix 
B(t)  is  a  (m  x  n)  distribution  matrix 
y_(t)  is  a  (p  x  1)  vector  of  system  outputs 
H(t)  is  a  (p  x  n)  measurement  matrix. 

It  will  be  assumed  that  measurement  of  the  output 
can  be  achieved  only  with  the  addition  of  noise,  so  the 
measurement  equation  is  defined  as 
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z(t)  =  y_(t)  +  v(t),  (2.3) 

where  v(t)  is  a  (p  x  1)  vector  representing  the  stochastic 
measurement  noise  and  z_(t)  is  a  (p  x  1)  vector  of  measured 
outputs . 

The  process  is  visualized  better  by  means  of  Fig. 2.1, 
where  the  double  lines  indicate  a  signal  flow  character- 
istic of  all  matrix  block  diagrams. 

It  is  assumed  for  simplicity  that  the  process  is  sta- 
tionary (time-invariant);   that  is,  all  the  elements  of 
the  A(t),  B(t)  and  H(t)  matrices  are  constants.   It  is 
further  assumed  that  the  process  is  to  be  sampled  at  dis- 
crete times  and  that  u(t)  and  w(t)  are  constant  vectors 
during  each  sampling  period.   The  general  solution  for 
Eq .  (2.1)  under  these  conditions  is  expressed  as 

x(t)  =  CAtx(t0)  +  CAtJc"ATB  [u  +  w]dT  .       (2.4) 

If  the  sampling  time  is  assumed  to  be  T  seconds,  Eq . 
(2.4)  may  be  rewritten  as 


x(T)  =  eATx(t0)  +  eAT  rTr"Ardr  bu  +  w]  .       (2.5) 

■J    n 


Defining 

3>(T)    =  rAT    ,  (2.6a) 
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A(t); 


Fig.  2.1  -  Block  diagram  for  a  continuous  stochastic 
process  with  time-varying  dynamics. 


i  ■■ 


v(k) 


z(k) 


Fig.  2.2  -  Block  diagram  for  a  discrete  stochastic 
process  with  time-invariant  dynamics. 
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and 


AT 


-An 


r(T)  ^  t:    /   - '  "  dT  B  , 
v  o 


(2.6b) 


where  $(T)  (n  x  n)  and  F  (T)  (n  x  m)  are  called  the  dis- 
crete state  transition  and  distribution  matrices,  respect- 
ively.  The  discretized   form  of  Eq.  (2.1)  becomes  in  the 
stationary  case 


x(k+l)  =  (j)(T)x(k)  +  r(T)|u(k)  +  w(k)l  .      (2.7) 


$(T)  and  T  (T)  may  be  obtained  by  application  of  the 
Laplace  transformation,  i.e., 


_  fi-l[ 


,-]] 


$(T)  =^"i!  (si  -  A)" 

L        J 


t-T  , 


(2.8a) 


Distinction  is  made  between  a  truly  discrete  system, and 
one  which  is  continuously  operating  but  driven  and  measured 
at  discrete  times.   The  linear  difference  equation  for  the 
truly  discrete  system  will  be 


x(k+l)  =  ADx(k)  +  BD  u(k)  +  w(k) 


The  terminology  "discrete",  "discretized"  is  used  without 
distinction  from  this  point  on. 


13 


r(T)  -^-1  liilS 


t  =T 


(2.8b) 


The  discretized,  linear  stochastic  process  with  time- 
invariant  coefficients  will  hereafter  be  represented  by 
the  model  (see  Fig.  2.2) 


x(k+l)  =  (jx(k)  +  r[u(k)  +  w(k) 


(2.9) 


with  measurement  equation 


z(k)  =  Hx(k)  +  v(k) 


(2.10) 


Eqs.  (2.9)  and  (2.10)  will  be  used  throughout  this 
thesis  as  the  model  of  the  complete  discrete  system  for 
which  reduced  sub-optimal  filters  are  implemented. 
Furthermore,  w(k)  and  v(k)  will  be  considered  to  be  vec- 
tors of  independent,  gaussian  random  processes  (white 
noise)  with  zero  means  and  fully  described  by  the  follow- 
ing statistics: 


Efw(k)l  =  0 


for  all  k  , 


(2.11) 


E  w(k)-wT(i)j  =  E[w(k)]  -E[wT(i)j  =0    k/  1,(2.12) 


r  T       I  ' 

Ejw(k)«w  (k)|  =  known  quantity  for  all  k  ,     (2.13) 


14 


E|v(k)(  =  0     for  all  k  , 


(2.14) 


E[v(k)-vT(i)j  =  E[v(k)]  -E[vT(i)j  =0    k  f   i,  (2.15) 


and 


E[v(k)-vT(k)i  =  R    for  all  k  , 


(2.16) 


where  R  is  a  (p  x  p)  matrix  of  covariance  of  measurement 
noise. 
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III.  IMPLEMENTATION  OF  A  REDUCED  DISCRETE, 
LINEAR  KALMAN  FILTER  BY  PARTITIONING 
THE  DYNAMICS  OF  THE  ORIGINAL  PROCESS 

This  chapter  presents  a  discrete  solution  to  the  prob- 
lem of  implementing  a  sub-optimal  filter  by  means  of  par- 
titioning the  original  discrete  system  (Fig.  2.2)  into  two 
reduced  models,  a  "primary  system"  model  which  will  include 
only  a  selected  subset  of  the  state  vector  components  whose 
estimation  is  of  primary  interest,  and  a  "secondary  sys- 
tem" which  will  be  formed  by  the  remaining  state  varibales. 
Then,  the  primary  plant  will  have  reduced  dimensions  and 
the  implementation  of  a  reduced  Kalman  filter  based  upon 
this  primary  model  will  be  possible. 

A.    DERIVATION  OF  THE  REDUCED  FILTER  EQUATIONS 

It  is  assumed  that  the  system  dynamics  are  represent- 
ed by  Eq .  (2.9)  and  that  the  statistics  of  the  forcing  and 
measurement  noise  are  fully  described  by  Eqs.  (2.11) 
through  (2.16).  For  simplicity  it  will  be  further  assumed 
that  there  is  no  deterministic  forcing  function  u(k).  The 
stochastic  process  may  then  be  expressed  as 

x(k+l)  =  (j)x(k)  +  Tw(k)  ,  (3.1) 

with  measurement 

z(k)  =  Hx(k)  +  v(k). 
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By  matrix  partitioning,  the  discrete  dynamics  of  the 
full  process  may  be  rewritten  as 


xL(k+l) 


x2(k+l) 


[*1  !  *cl 


I 


c2  ,  *2 


xx(k) 


x2(k) 


+ 


n 


^  I  o  [wx(k)] 


o    '  rJ  w0(k) 


(3.2) 


I  Zll~Z        J 


With  the  full  state  vector  being  an  (n  x  1)  matrix, 
then : 

X-,  represents  an  (n-,  x  1)  state  vector  which  includes 
all  those  states  of  particular  interest  in  the  estimation 
process ; 

Xy    represents  an  (n~  x  1)  vector  of  the  remaining 
state  vector  components  (n-,  +  n~  =  n); 

(j)-,  is  an  (n-,  x  n-,  )  matrix,  representing  the  dynamics 
of  X-j^ 

$2  is  an  (n2  x  n^)  matrix,  representing  the  dynamics 
of  Xo  ; 

$  -,  and  I  2  are  (ni  x  no)  an(^  ^n2  x  nl  ^  matrices  re- 
spectively, indicating  the  coupling  between  the  primary 
and  secondary  systems; 

w-,  is  a  (m-.  x  1)  vector  of  forcing  (white  noise) 
functions ; 

Wo  is  a  (tt)2  x  1)  vector  of  forcing  (white  noise) 
functions  with  m-,  +  m2  =  m;  and 

I~^  and  T\    are  (n.  x  m-,  )  ,  (n~  x  m~)  distribution  ma- 
trices respectively. 
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The  observation  matrix  H  is  partitioned  as  follows, 


H  =   H,  I  H2  j  . 


The  measurement  equation  may  therefore  be  rewritten  as 

z(k)  =  |  HL  ;  H2  j  [  X]_(k)  ]  +  v(k)  .  (3.3) 

[i27kTj 

It  is  assumed  that  all  the  measured  outputs  are  in- 
cluded in  the  primary  system,  which  makes  H2  =  0.   This 
partitioning  may  not  always  be  possible,  and  even  if  it 
is,  it  may  require  special  rearrangement  of  the  z   vector. 

In  Eq .  (3.3)  z^(k)  is  the  original  (p  x  1)  measured 
vector  and  H-,  is  the  (p  x  n-,  )  measurement  matrix  for  the 
primary  system.   It  may  be  seen  from  this  special  parti- 
tioning of  the  H  matrix  that  the  primary  system  vector 
X-.  will  not  only  be  formed  by  those  states  of  primary  in- 
terest in  the  estimation  process,  but  will  also  include 
those  states  whose  estimation  may  not  be  needed,  but 
whose  outputs  are  being  measured.   This  will  establish 
the  condition  that  n-,*p. 

As  a  result  of  the  preceeding  discussion,  Eq .  (3.3) 
becomes 

z(k)  =  H1x1(k)  +  v(k)  ,  (3.4) 
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which  represents  the  measurement  equation  for  the  reduced 
system. 

The  general  system  may  be  described  as  consisting  of 
a  primary  system. 

x^k+1)  -  ^x-^k)  +  $clx2(k)  +  rjw^k)  ,     (3.5a) 

and  a  secondary  system  described  by, 

x2(k+l)  =  $c22ii(k)  +  $2x2(k)  +  T2w2(k)  .     (3.5b) 

Fig.  (3.1)  shows  a  complete  representation  of  the  parti- 
tioned system. 

In  order  to  reduce  computational  requirements  in  the 
estimation  of  the  preselected  states,  the  design  of  the 
filter  is  based  on  the  dynamics  of  the  primary  system, 
which  can  be  represented  as  a  simplified  model  of  the  re- 
duced process  by 

xL(k+l)  =  ^x^k)  +  Fjw^k)  ,  1  (3.6) 


1 

A  bar  over  a  quantity  indicates  a  reduced-model  quan' 
tity. 
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with  measurement 


z(k)  =  H^Ck)  +  v(k)  . 


In  the  simplified  model  of  the  primary  system  the  term  in 
volving  $  i2£?(k)  has  been  neglected  due  to  the  fact  that 
the  secondary  state  vector  x~(k)  is  not  considered  at  all 
in  the  implementation  of  the  reduced  filter.   The  conse- 
quences of  using  Eq .  (3.6)  instead  of  Eq .  (2.5a),  will  be 
analyzed  in  Section  B  of  this  chapter. 

If  the  Kalman  filter  algorithm  is  applied  using  the 
dynamics  of  the  full  plant,  the  matrix  equation  for  the 
gain  G(k),  the  covariance  of  estimation  error  P(k/k)  and 
the  prediction  covariance  P(k+l/k)  are  given  in  the  dis- 
crete case  by 

G(k)  =  P(k/k-l)HT   H  P(k/k-l)HT  +  R  "   ,     (3.7a) 


P(k/k)  =   I  -  G(k)H  P(k/k-l)  , 


(3.7b) 


and 


P(k+l/k)  =  $  P(k/k)$T  +  Q  , 


(3.7c) 


where 


Q  =  TE[w(k)wT(k) 


r 


G(k)  is  an  (n  x  p)  matrix,  P(k/k),  P(k/k-l)  and  P(k+l/k) 
are  (n  x  n)matrices. 
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The  estimation  equation  is  given  by 


x(k/k)  =  $x(k-l/k-l)  +  G(k)fz(k)  -  H^xCk-l/k-l)1 


(3.8) 


where  x(k/k)  is  an  (n  x  1)  vector  of  estimated  states. 
By  induction,  the  reduced  sub-optimal  expressions 
using  the  dynamics  of  Eq .  (3.6)  are 


-1 


GL(k)  =  P1(k/k-l)H^[H1P1(k/k-l)H^  +  It]    ,     (3.9a) 


P1(k/k)  =  'i  -  G1(k)H1.P1(k/k-l) 


(3.9b) 


and  similarly 


P1(k+l/k)  =  $1P1(k/k)$J  +  Ql    , 


(3.9c) 


where 


Qi  -  ^E 


w1(k)w^(k)]r[ 


l(k/k)  (n-,  x  n,  )  represents  the  sub-optimal  covariance 
of  estimation  error  and  P-,  (k+l/k)  (n-.  x  n,  )  is  the  sub- 
optimal  prediction  covariance  of  estimation  error. 
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The  sub-optimal  estimation  equation  is 


^(k/k)  =  (jjjc^k-l/k-l)  +  G]_(k)[z(k)  -  H1$1x1(k-l/k-l)]  . 


(3.10) 

Considering  the  problem  of  initialization  for  the 
discrete  reduced  covariance  equations,  it  is  seen  as  a 
good  approach  to  use  for  P-.(0/-l),  the  partitioned  por- 
tion of  P(0/-1)  corresponding  to  P,-|(0/-1),  or  considering 
that 


P(0/-1)   = 


Pn(o/-i)      p12(o/-i) 


p21(o/-i)      p22(o/-i) 


then 


P1(0/-1)  =  p11(o/-i)   . 


For  the  filter  initial  states, 


where 


x^O/-!)  =  x^O/-!)  , 


M 


>- 


x(0/-l)   = 


I2"(o7-~i)~ 
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Fig.  3.2  shows  a  complete  block  diagram  of  the  reduced 
plant  and  filter. 

B.    REDUCED  FILTER  DEGRADATION 

Because  the  dynamics  of  a  simplified  system  repre- 
sented by 

x-^k+1)  =  (Jl^Ck)  +  IjWjXk)  (3.6) 

were  used  for  the  implementation  of  the  reduced  filter, 
instead  of  the  dynamics  of  the  primary  system 

xx(k+l)  =  ^x-^k)  +  ^clx2(k)   +   rjw-^k)  ,      (3.5a) 

the  estimation  of  x-,  is  caused  to  be  sub-optimal  and  in 
error.   This  error  may  be  considered  the  performance  de- 
gradation of  the  sub-optimal  filter. 

Considering  the  reduced  filter,  its  estimation  error 
covariance  is  given  by 


V(k)  =  El  x^k)  -  Xi(k/k) 

i  L   J.  J_  j 


(3.11) 


where 


xx(k)  =  ^x^k-l)  +  <£clx2(k-l)  +  Fjw^k-1)  ,  (3.5a) 

and  by  proper  use  of  Eq .  (3.6)  it  may  be  rewritten  as 

xx(k)  =  xx(k)  +  $clx2(k-l)  .  (3.12) 
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In  the  same  way, 


x^k/k)  =  $1|1(k-l/k-l)+G1(k)[z(k)  -  H^x^k-l/k-1)]  , 

(3.10) 
which  by  substitution  of  Eq .  (3.4)  for  z(k),  gives 
xx(k/k)  =  $1x1(k-l/k-l)+  G1(k)[H1x1(k)+v(k)-H1$1|1(k-lA-D 


Further  substitution  of  Eq .  (3.12)  for  x-,(k:)  gives 

xx(k/k)  =  (^(k-l/k-l)  +  G1(k)H1|x1(k)  +  $clx2(k-l) 


-  (^(k-l/k-1) 


+  Gx(k)v(k)  . 


(3.13) 


Using  Eqs.  (3.12)  and  (3.13),  the  experimental  estimation 
error  covariance  V(k)  (Eq.  3.11)  may  be  rewritten,  after 
proper  rearrangement,  as 

r  .- 

V(k)  =  E  l(Ix  -  G1(k)H]_)(x1(k)  -  (Jjc^k-l/k-l) 

+  (l1-G1(k)H1)$clx2(k-l)  +  G1(k)v(k)j  |. 

(3.14) 

By  application  of  matrix  algebra  and  elemental  algebra  of 
expectations,  Eq .  (3.14)  may  be  expanded  and  after  some 
manipulation  gives 
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Vi 


(k)  =  [l1-G1(k)H1]P1(k/k-l)+[l1-G1(k)H1]   ^[x^k-l)]  $^ 


r  2 


+E||1(k/k-l)x2(k-l)]  ^el+^c^t^'1^^"1^  [Ii"Gi(k)H 


(3.15) 


where 


/\ 


e^k/k-1)  =  xx(k)  -  x^k-l/k-1)  . 

The  first  term  on  the  right  hand  side  of  Eq .  (3.15) 
can  be  easily  recognized  as  P-,(k/k),  the  sub-optimal 
theoretical  variance  of  estimation  error  based  on  the  re- 
duced model.   It  is  intuitively  expected  that  P-.(k/k)  > 
P-|-,(k/k).   The  remaining  terms  of  Eq .  (3.15)  represent  the 
further  degradation  of  the  reduced  filter  and  are  seen  to 
be  functions  of  that  partitioned  portion  of  ^  which  couples 
Xy    to  the  primary  system,  i.e.,  $„]•   This  suggests  that  as 
the  coupling  between  Xo  an<^  tn-e  primary  system  ($  -,  )  de- 
creases, the  quantity  P-.  (k/k)  more  accurately  represents 
the  actual  estimation  error  covariance. 
C.    EXAMPLE  1 

The  process  used  for  this  example  was  a  fourth-order 
model,  with  continuous  linear  stationary  dynamics  subject- 
ed to  a  forcing  function  having  both  a  deterministic  and  a 
random  component,  with  the  latter  having  zero  mean  and 
known  variance.   It  was  assumed  that  only  one  state  vari- 
able was  measured,  subject  to  the  presence  of  gaussian 
measurement  noise  with  zero  mean  and  known  variance.   It 
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was  further  assumed  that  both  random  sequences  have  sta- 
tistics as  described  by  Eqs.  (2.11)  through  (2.16).  The 
reduced  filter  was  implemented  in  order  to  estimate  only 
two  states  (x-,,X2). 

The  fourth-order  continuous  system  is  described  by 


x(t)  = 


0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

400 

-320 

-118 

-19 

x(t)  + 


0 

0 

0 

400 


u(t)+w(t) 


A  deterministic  function  u(t)  =0.5  was  used  in  the  ex- 
ample. The  random  forcing  function  w(t)  was  assumed  to 
have  zero  mean  and  variance  equal  to  0.01  square  units, 
i.e. ,  w(t)  :  N(0.0,  0.01). 

The  measurement  equation  is 

z(t)  =  [l   0   0   o]x(t)  +  v(t)  , 

with  the  measurement  noise  v(t),  considered  to  have  zero 
mean  and  variance  equal  to  0.25  square  units,  i.e., 
v(t)  :  N(0.0,  0.25). 

Before  applying  the  reduced-filter  technique  dis- 
cussed in  Section  III,  A,  the  continuous  process  was  dis- 
cretized  using  the  ideas  presented  in  Chapter  II,  so  the 
discretized  system  could  be  described  as 

x(k+l)  =  $x(k)  +  rru(k)  +  w(k) 
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with  measurement 


z(k)  = 


10   0   0 


x(k)  +  v(k)  , 


where 


$  -  i  AT 


(2.6a) 


was  computed  using  the  series  expansion 


2    i.3 


$  =  I  +  AT  +  yr(AT)  +  Jr(AT)  + 


2.' 
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and 


r  =  eAT  rTfc-Ardr  B  , 

J   o 


(2.6b) 


was  computed  by  means  of 


I  =  (I  +  AT  +  yr(AT)  + 


?        3 
)(IT  -  A  i-,-  +  A z  ^-  -...)B, 


A  Monte  Carlo  simulation   was  performed  using  200 
replications  of  the  trajectory,  each  trajectory  being 
formed  by  50  samples  with  sampling  time  T  =  0.05  seconds. 


An  IBM  System/360-67  computer  was  used  for  solving  all 
the  examples  presented  in  this  thesis.   A  copy  of  the 
main  computer  program  for  this  example  is  presented  in 
Appendix  B. 
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This  simulation  was  carried  out  in  order  to  produce  approx- 
imate ensemble  values  for  the  experimental  mean  and  experi- 
mental variance  of  estimation  error  in  x-,  and  x^ ,  as  a  func- 
tion of  time  from  t  =  0,  for  both  the  complete,  optimal  fil- 
ter and  the  reduced  filter. 

For  the  Monte  Carlo  simulation  each  state  was  assumed 
to  have  uncorrelated  gaussian  initial  conditions  with  zero 
mean  and  the  following  variances: 


(T2x 

J'2x 


-,  =  0.5  square  units 
2  =  2.0  square  units 


0    Xo  =  2.0  square  units 


cf2 


x,  =  4.0  square  units 


In  order  to  obtain  an  unbiased  estimate  from  the 
Kalman  filter,  the  initial  values  assigned  to  the  filter 
state  variables  were 


£(o/-i)  =  e[x(o)] 


or 


x(0/-l)   = 


0 
0 
0 
0 


30 


md  the  initial  value  of  P(0/-1)  was 


P(0/-1)  =  E;  ;x(0)  -  x(0/-l)j  [x(0)  -  x(0/-l)] 

r. 


=  E{ 


x(0)  -  Ejx(0)   jx(0) 


-  E 


x(0) 


=  cov  x(0) |  , 


or 


P(0/-1)  = 


0.5 


2.0 


O 


o 


2.0 


L 


4.0 


The  reduced  filter  was  implemented  to  estimate  only 
two  states  (x-,  and  X2),  and  for  this  purpose  partitioning 
was  applied  as  follows.   The  $  and  V    matrices  (from 
Table  III-l) 


$(0.05)  = 


0.9990   0.0449  '  0.0012    0.0000 

l 

-0.0066    0.9946  I  0.0480   0.0009 


-0.3644   -0.2981  I  0.8871    0.0307 

I 
-12.2713  -10.1815  1-3.9181    0.3043 


rT(0.05)  =  fo.0001   0.0066  !  0.3644    12.2713], 
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were  partitioned  as  shown  to  obtain 


4i  " 


0.9990    0.0449 
-0.0066    0.9946 


and 


r,  = 


0.0001 
0.0066 


Correspondingly, 


P1(0/-1)  = 


0.5    0 
0    2.0 


x^O/-!)  = 


L°J 


H,  =  i  1   0 


and 


Qi  = 


0.000   0.000  " 


0.000   0.000 


.   (See  Table  III-l) 
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The  reduced  filter  was  then  implemented  using  the  set 

of  Eqs.  (3.9a,b,c)  repeated  here  for  convenience, 

0000.0    n ;  •  •  ■ 

G1(k)  =  P1(k/k-l)HJ['H1P1(k/k-l)H1  +  R     ,    (3.9a) 


i 


, 


. 


P1(k/k)  =  [I  -  G1(k)H1jP1(k/k-l)  , 


(3.9b) 


and 


P1(k+l/k) 


=  $]_P1(k/k)$1  +  Q 


(3.9c) 


Attention  should  be  called  to  the  fact  that  the  pre- 
sence of  a  deterministic  forcing  function  in  this  problem 
changes  slightly  the  form  of  £Qjv  (3.10),  which  can  be  re- 
written  as 


x-^k/k)  =  |1f1(k-l/k-l)+  Eu  (k-l)+  G]_(k) 


(k) 


Tt  Cl 


-H1:$1x1(k-l/k-l)+  I^uj^k-l) 


to  reflect  the  effect  .of  u(k)  on  the  prediction  step. 

The  outputs  -of  the  computer  simulation  are  shown  in 
Tables  III-2  through  III-5,  where  numerical  results  versus 
time  are  tabulated  for  the  gains,  theoretical  variance, 
experimental  variance  and  experimental  means  of  estimation 
error  for  both  the  complete  and  reduced  filters.   Figs. 
(3.3)  through  (3.6)  show  graphical  representations  com- 
paring complete  and  reduced  filters  performance. 
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Table  III-l 

PHI  MATRIX 

0.9999 

0.0499 

0.0012 

0.0000 

-0.0066 

0.9946 

0.0480 

0.0009 

-0.3644 

-0.2981 

0.8871 

0.0307 

-12.2713 

-10.1815 

-3.9181 

0.3043 

GAMMA  MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 

Q  MATRIX 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0008 

0.0000 

0.0000 

0.0013 

0.0447 

0.0000 

0.0008 

0.0447 

1.5059 

VARIANCE  OF  FORCING  NOISE  •  0.01 


VARIANCE  OF  MEASUREMENT  NOISE  -  0.25 


INITIALIZATION  -  P(0/-1) 


0.5000 

0.0 

0.0 

0.0 

0.0 

2.0000 

0.0 

0.0 

0.0 

0.0 

2.0000 

0.0 

0.0 

0.0 

0.0 

4.0000 

INITIALIZATION  -  X(0/-1) 

0.0 
0.0 
0.0 
0.0 


T  =  0.05  seconds. 
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Table  III-2 
GAINS  FOR  FULL  AND  REDUCED  FILTERS 


K 

G<  1,1) 

G(2,ll 

G(3, 1) 

G(4,l) 

GR< 1,11 

GR(2,1I 

1 

0.6667 

0.0 

0.0 

0.0 

0.6667 

0.0 

2 

0.4071 

0.2333 

-0.2095 

-7.2845 

0.4071 

0.2330 

3 

0.3097 

0.4075 

-0.6537 

-9.4633 

0.3102 

0.4259 

4 

0.2677 

0.5088 

-1.1719 

-9.0573 

0.2708 

0.  5779 

5 

0.2477 

0.5314 

-1.6080 

-7.0288 

0. 2569 

0.6814 

6 

0.2349 

0.4P01 

-1.8640 

-4.4015 

0.2540 

0.7347 

7 

0.2227 

0.4108 

-1.9266 

-1.9563 

0. 2542 

0.7445 

8 

0.2087 

0.3175 

-1.8378 

-0.0721 

0.2537 

0.7222 

q 

0. 1928 

0.2268 

-1.6557 

1.1851 

0.2512 

0.6798 

10 

0.175P 

0.1474 

-1.4305 

1.9077 

0. 2464 

0.6274 

11 

0.  1584 

0.0825 

-1.1970 

2.2348 

0.2399 

0.5720 

12 

0.1412 

0.0322 

-0.9760 

2.2961 

0.2321 

0.5177 

13 

0.1247 

-0.0051 

-0.7780 

2.1931 

0.2237 

0.4669 

14 

0.  1093 

-0.0315 

-0.6072 

1.9984 

0.2151. 

0.4205 

15 

0.0950 

-0.0490 

-0.4637 

1.7610 

0.2064 

0.3789 

16 

0.0819 

-0.0596 

-0.3457 

1.5120 

0.1979 

0.3418 

17 

0.0701 

-0.0647 

-0.2505 

1.2705 

0.1897 

0.3090 

18 

0.0596 

-0.0660 

-0.1752 

1.0474 

0.1819 

0.2800 

19 

0.0502 

-0.0643 

-0.1167 

0.8480 

0.1745 

0.2544 

20 

0.0420 

-0.0606 

-0.0723 

0.6746 

0. 1675 

0.2317 

21 

0.0348 

-0.0556 

-0.0395 

0.5269 

0. 1609 

0.2115 

22 

0.0286 

-0.0499 

-0.0161 

0.4036 

0.1547 

0.1937 

23 

0.0233 

-0.0438 

-0.0001 

0.3026 

0.  1488 

0.1777 

24 

0.0188 

-0.0378 

0.0100 

0.2215 

0.1433 

0.1635 

25 

0.0151 

-0.0320 

0.0157 

0.1576 

0.1382 

0.1508 

26 

0.0120 

-0.0266 

0.0181 

0.1084 

0. 1333 

0.1393 

27 

0.0095 

-0.0217 

0.0181 

0.0714 

0.1287 

0.1290 

28 

0.0075 

-0.0174 

0.0166 

0.0444 

0. 1243 

0.1197 

29 

0.0059 

-0.0136 

0.0141 

0.0254 

0.1203 

0.  1112 

30 

0.0047 

-0.0105 

0.0111 

0.0127 

0.1164 

0.1036 

31 

0.0037 

-0.0078 

0.0080 

0.0046 

0.1127 

0.0966 

32 

0.0031 

-0.0057 

0.0050 

-0.0000 

0.1092 

0.0902 

33 

0.0026 

-0.0041 

0.0023 

-0.0022 

0.1059 

0.0844 

34 

0.0022 

-0.0028 

-0.0001 

-0.0028 

0.1028 

0. 0790 

35 

0.0020 

-0.0018 

-0.0022 

-0.0024 

0.0998 

0.0741 

36 

0.0018 

-0.0011 

-0.0038 

-0.0015 

0.0970 

0.0696 

37 

0.0018 

-0.0006 

-D.0050 

-0.0004 

0.0943 

0.0654 

38 

0.0017 

-0.0003 

-0.0060 

0.0007 

0.0917 

0.0616 

39 

0.0017 

-0.0001 

-0.0066 

0.0017 

0. 0893 

0.0580 

40 

0.0017 

0.0000 

-0.0070 

0.0024 

0.0869 

0.0547 

41 

0.0017 

0.0001 

-0.0073 

0.0030 

0.0847 

0.0517 

42 

0.0017 

0.0001 

-0.0074 

0.0033 

0. 0825 

0.0488 

43 

0.0017 

0.0001 

-0.0074 

0.0034 

0.0804 

0.0461 

44 

0.0017 

0.0001 

-0.0073 

0.0033 

0.0785 

0.0437 

45 

0.0017 

0.0000 

-0.0072 

0.0031 

0.0765 

0.0413 

46 

0.0017 

0.0000 

-0.0071 

0.0028 

0.0747 

0.0392 

47 

0.0017 

-0.0000 

-0.0070 

0.0025 

0.0730 

0.0371 

48 

0.0017 

-0.0000 

-0.0068 

0.0022 

0.0713 

0.0352 

49 

0.0017 

-0.0001 

-0.0067 

0.0018 

0.0696 

0.0334 

50 

0.0017 

-0.0001 

-0.0066 

0.0015 

0.0680 

0.0318 

G  (4x1) 

represents 

the  matrix  of  gains 

for  the 

optimal 

filter. 

GR(2xl) 

represents 

the  matrix  of  gains 

for  the 

reduced 

filter. 
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Table  III-3 
THEORETICAL  VARIANCE  OF  ESTIMATION  ERROR 


K 

PK<1,1) 

PK(2t2l 

PM3,3) 

PK(4,4) 

PR(1,1> 

PR(2t2l 

1 

0.1667 

2.0000 

2.0000 

4.0000 

0.1667 

2.0000 

2 

0.1018 

1.9603 

1.7605 

242.6279 

0.1018 

1.9558 

3 

0.0774 

1.8002 

2.3189 

346.3362 

0.0776 

1.8684 

4 

0.0669 

1.5182 

3.9792 

261.1479 

0.0677 

1.7325 

5 

0.0619 

1.1771 

5.7036 

142.2309 

0.0642 

1.5559 

6 

0.0587 

0.8478 

6.6569 

60.6997 

0.0635 

1.3561 

7 

0.0557 

0.5745 

6.6825 

21.3209 

0.0635 

1.1534 

8 

0.0522 

0.3709 

6.0419 

8.2223 

0.0634 

0.9640 

9 

0.0482 

0.2307 

5.0782 

6.9855 

0.0628 

0.7970 

10 

0.0439 

0.1400 

4.0514 

9.3412 

0.0616 

0.6557 

11 

0.0396 

0.0847 

3.1097 

11.6978 

0.0600 

0.5391 

12 

0.0353 

0.0529 

2.3160 

12.9672 

0.0580 

0.4442 

13 

0.0312 

0.0361 

1.6824 

13.1242 

0.0559 

0.3676 

14 

0.0273 

0.0282 

1.1960 

12.4778 

0.0538 

0.3058 

15 

0.0237 

0.0255 

0.8337 

11.3650 

0.0516 

0.2559 

16 

0.0205 

0.0254 

0.5710 

10.0516 

0.0495 

0.2155 

17 

0.0175 

0.0264 

0.3855 

8.7171 

0.0474 

0.1827 

18 

0.0149 

0.0275 

0.2581 

7.4696 

0.0455 

0.1557 

19 

0.0126 

0.0283 

0.1736 

6.3652 

0.0436 

0.1336 

20 

0.0105 

0.0286 

0.1199 

5.4260 

0.0419 

0.1152 

21 

0.0087 

0.0282 

0.0877 

4.6526 

0.0402 

0.0999 

22 

0.0072 

0.0273 

0.0701 

4.0336 

0.0387 

0.0870 

23 

0.0058 

0.0259 

0.0619 

3.5512 

0.0372 

0.0762 

24 

0.0047 

0.0241 

0.0595 

3.1854 

0.0358 

0.0670 

25 

0.0038 

0.0221 

0.0601 

2.9159 

0.0345 

0.0591 

26 

0.0030 

0.0199 

0.0619 

2.7237 

0.0333 

0.0524 

27 

0.0024 

0.0177 

0.0639 

2.5918 

0.0322 

0.0466 

28 

0.0019 

0.0155 

0.0654 

2.5057 

0.0311 

0.0416 

29 

0.0015 

0.0133 

0.0659 

2.4529 

0.0301 

0.0373 

30 

0.0012 

0.0114 

0.0654 

2.4236 

0.0291 

0.0335 

31 

0.0009 

0.0096 

0.0640 

2.4099 

0.0282 

0.030? 

32 

0. 00C8 

0.0081 

0.0618 

2.4059 

0.0273 

0.027? 

33 

0.0006 

0.0067 

0.0590 

2.4071 

0.0265 

0.0247 

34 

0.0006 

0.0056 

0.0557 

2.4105 

0.0257 

0.0224 

35 

0.0005 

0.0046 

0.0523 

2.4141 

0.0250 

0.0204 

36 

0.0005 

0.0039 

0.0489 

2.4167 

0.0243 

0.0186 

37 

0.0004 

0.0033 

0.0455 

2.4177 

0.0236 

0.0170 

38 

0.0004 

0.0028 

0.0424 

2.4170 

0.0229 

0.0155 

39 

0.0004 

0.0024 

0.0396 

2.4146 

0.0223 

0.0142 

40 

0.0004 

0.0022 

0.0371 

2.4109 

0.0217 

0.0131 

41 

0.0004 

0.0020 

0.0349 

2.4062 

0.0212 

0.0120 

42 

0.0004 

0.0018 

0.0330 

2.4008 

0.0206 

0.0111 

43 

0.0004 

0.0017 

0.0315 

2.3951 

0.0201 

0.0102 

44 

0.0004 

0.0017 

0.0303 

2.3894 

0.0196 

0.0094 

45 

0.0004 

0.0017 

0.0293 

2.3839 

0.0191 

0.0037 

46 

0.0004 

0.0016 

0.0285 

2.3787 

0.0187 

0.0081 

47 

0.0004 

0.0016 

0.0279 

2.3741 

0.0182 

0.0075 

48 

0.0004 

0.0016 

0.0275 

2.3700 

0.0178 

0.0070 

49 

0.0004 

0.0016 

0.0272 

2.3664 

0.0174 

0.0065 

50 

0.0004 

0.0016 

0.0270 

2.3635 

0.0170 

0.0060 

PK  represents  the  diagonal  elements  of  the  theoretical 

covariance  matrix  for  the  optimal  filter. 

PR  represents  the  diagonal  elements  of  the  theoretical 

covariance  matrix  for  the  reduced  filter. 
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Table    III-4 
EXPERIMENTAL    VARIANCE    OF    ESTIMATION    ERROR 

K     VK(1,1)    VK<2,2)    VM3.3)    VK<4,4)    VRdtll    VP(2,2» 


1 
2 
3 

4 
5 
6 

7 
8 

Q 

10 

11 

12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 


0.1510 
0.1022 
0.0807 
0.0667 
0.0677 
0.0529 
0.0518 
0.0487 
0.0464 
0.0457 
0.0420 
0.0391 
0.0368 
0.0297 
0.0264 
0.0226 
0.0188 
0.0153 
0.0131 
0.0112 
0.0094 
0.0074 
0.0061 
0.0047 
0.0039 
0.0032 
0.0025 
0.0020 
0.0016 
0.0013 
0.0010 
0.0008 
0.0007 
0.0006 
0.0005 
0.0005 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 
0.0004 


1.9066 

1.8632 

1.7518 

1.5112 

1.1743 

0.8372 

0.5928 

0.3939 

0.2513 

0.1612 

0.1002 

0.0617 

0.0401 

0.0300 

0.0256 

0.0247 

0.0260 

0.0261 

0.0271 

0.026  9 

0.0274 

0.0252 

0.0240 

0.0216 

0.0202 

0.0185 

0.0166 

0.0149 

0.0130 

0.0114 

0.0098 

0.0083 

0.0071 

0.0060 

0.0050 

0.0042 

0.0035 

0.0029 

0.0025 

0.0022 

0.0019" 

0.0018 

0.0017 

0.0016 

0.0015 

0.0014 

0.0014 

0.0014 

0.0015 

0.0016 


1.9733 

1.7310 
2.2373 
3.8145 
5.6299 
6.0810 
6.5273 
6.0208 
5.1600 
4.4359 
3.5490 
2.7683 
2.0579 
1.3831 
0.9971 
0.6811 
0.4289 
0.2742 
0.1830 
0.1292 
0.0892 
0.0707 
0.0651 
0.0634 
0.0596 
0.0569 
0.0547 
0.0546 
0.0550 
0.0581 
0.0584 
0.0567 
0.0539 
0.0490 
0.0435 
0.0420 
0.0437 
0.0446 
0.0421 
0.0366 
0.0300 
0.0256 
0.0241 
0.0261 
0.0300 
0.0317 
0.0307 
0.0299 
0.0281 
0.0229 


4. 
225. 
326. 
257. 
138. 
60. 
23. 

8. 

7. 

9. 
12. 
14. 
14. 
14. 
11. 
12. 

9. 

7. 

6. 

5. 

5. 

4. 

3. 

2. 

3. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

2. 

1. 

1. 

2. 

2. 

2. 

2. 


0880 
4341 
4941 
3462 
0874 
2167 
2761 
7271 
2510 
6855 
4134 
0120 
2776 
0637 
9529 
3539 
7735 
4905 
0835 
5882 
2477 
1038 
5547 
8311 
2233 
5204 
7723 
8033 
6636 
0651 
4213 
7746 
5252 
1855 
5705 
2536 
5342 
3070 
0512 
4148 
2901 
5146 
2052 
3978 
9640 
9171 
2013 
3214 
2047 
8952 


0.1486 
0. 1008 
0.0794 
0.0652 
0.066? 
0.0511 
0.0504 
0.0492 
0.0497 
0.0542 
0. 0570 
0.0559 
0.0615 
0.0592 
0.0587 
0.0593 
0.0585 
0.0616 
0.0639 
0.0661 
0.0694 
0.0661 
0.0661 
0.0649 
0.0640 
0.0593 
0.0555 
0.0566 
0.0542 
0.0510 
0.0474 
0.0459 
0.0435 
0.0391 
0.0378 
0.0363 
0.0349 
0.0350 
0.0378 
0.0393 
0.0403 
0.0405 
0.0411 
0.0429 
0.0439 
0.0476 
0.0494 
0.0500 
0.0532 
0.0570 


1.9041 
1.8617 
1.7496 
1.5140 
1.2011 
0.9078 
0.7632 
0.7093 
0.7732 
0.9210 
1.0924 
1.1048 
1.1834 
1.2549 
1.2175 
1.1740 
1.0867 
1.0225 
0.9195 
0.8078 
0.7044 
0.5994 
0.4966 
0.4195 
0.3411 
0.2732 
0.2257 
0.1949 
0.1729 
0.1560 
0.1489 
0.1498 
0.1523 
0.1563 
0.1674 
0.1752 
0.1817 
0.1891 
0.1947 
0.1985 
0.2024 
0.2035 
0.2022 
0.2016 
0.19O0 
0.1957 
0.1914 
0.1835 
0.1772 
0.1698 


VK  represents  the  diagonal  elements  of  the  experimental 
covariance  matrix  for  the  optimal  process. 

VR  represents  the  diagonal  elements  of  the  experimental 
covariance  matrix  for  the  reduced  process. 


37 


Table  III-5 
EXPERIMENTAL  MEANS  OF  ESTIMATION  ERROR 


K 

MEANE1 

MEANE2 

MEANE3 

MEANE4 

MEANE1R 

MEANE2R 

1 

0.0496 

0.0500 

-0.0350 

-0.1543 

0.0496 

0.0500 

2 

0.0372 

0.0392 

-0.0572 

-0.6276 

0.0372 

0.0392 

3 

0.0356 

0.0307 

-0.0903 

-0.7992 

0.0358 

0.0429 

4 

0.0380 

0.0271 

-0.1340 

-0.7676 

0.0391 

0.0662 

5 

0.0395 

0.0203 

-0.1637 

-0.3713 

0.0425 

0.0973 

6 

0.0333 

-0.0030 

-0.1219 

-0.0534 

0.0387 

0.1099 

7 

0.0299 

-0.0146 

-0.0909 

0.2345 

0.0387 

0.1346 

8 

0.0044 

-0.0565 

0.1352 

0.1478 

0.0122 

0.0810 

9 

-0.0017 

-0.0536 

0.1711 

0.0767 

0.0089 

0.1003 

10 

0.0049 

-0.0373 

0.1013 

0.1589 

0.0229 

0.1579 

11 

0.0142 

-0.0265 

0.0168 

-0.0124 

0.0414 

0.2121 

12 

0.0021 

-0.0281 

0.0951 

-0.0444 

0.0257 

0.1751 

13 

-0.0028 

-0.0232 

0.1136 

-0.1564 

0.0209 

0.1623 

14 

-0.0072 

-0.0168 

0.1219 

-0.3115 

0.0155 

0.1447 

15 

-0.0086 

-0.0108 

0.1078 

-0.4010 

0.0148 

0.1327 

16 

-0.0058 

-0.0083 

0.0729 

-0.3752 

0.0232 

0.1322 

17 

-0.0028 

-0.0080 

0.0511 

0.0015 

0.0320 

0.1267 

18 

0.0020 

-0.0112 

0.0313 

-0.0783 

0.0461 

0.1255 

19 

0.0055 

-0.0150 

0.0153 

-0.0996 

0.0572 

0.1156 

20 

0.0044 

-0.0137 

0.0109 

-0.1052 

0.0511 

0.0790 

21 

0.0040 

-0.0136 

0.0142 

0.2029 

0. 0476 

0.0468 

22 

0.0014 

-0.0094 

0.0236 

0.1062 

0.0321 

0.0012 

23 

0.0022 

-0.0103 

0.0273 

0.0610 

0.0343 

-0.0198 

24 

0.0017 

-0.0089 

0.0262 

-0.0824 

0.0284 

-0.0492 

25 

0.0011 

-0.0072 

0.0235 

-0.0285 

0.0199 

-0.0795 

26 

0.0008 

-0.0063 

0.0215 

-0.0504 

0.0143 

-0.1043 

27 

0.0007 

-0.0056 

0.0207 

-0.0033 

0.0100 

-0.1258 

28 

0.0006 

-0.0048 

0.0263 

0.1783 

0.0050 

-0.1455 

29 

0.0001 

-0.0027 

0.0329 

0.1111 

-0.0078 

-0.1700 

30 

0.0002 

-0.0014 

0.0292 

-0.2192 

-0.0114 

-0.1832 

31 

0.0003 

-0.0006 

0.0215 

-0.1151 

-0.0128 

-0.1929 

32 

0.0004 

0.0002 

0.0170 

-0.0758 

-0.0172 

-0.2037 

33 

0.0005 

0.0008 

0.0080 

-0.2516 

-0.0234 

-0.2143 

34 

0.0004 

0.0011 

0.0012 

-0.0467 

-0.0344 

-0.2271 

35 

0.0004 

0.0012 

0.0022 

0.0639 

-0.0455 

-0.2379 

36 

0.0005 

0.0014 

0.0016 

-0.0675 

-0.0514 

-0.2430 

37 

0.0006 

0.0014 

-0.0011 

-0.0404 

-0.0574 

-0.2467 

38 

0.0006 

0.0012 

-0.0069 

-0.1685 

-0.0654 

-0.2506 

39 

0.0006 

0.0007 

-0.0093 

0.0404 

-0.0711 

-0.2520 

40 

0.0006 

0.0003 

-0.0057 

0.0877 

-0.0785 

-0.2533 

41 

0.0007 

0.0002 

-0.0006 

0.1207 

-0.0795 

-0.2494 

42 

0.0007 

0.0003 

0.0027 

0.0151 

-0.0876 

-0.2489 

43 

0.0008 

0.0004 

0.0038 

0.0374 

-0.0882 

-0.2432 

44 

0.0008 

0.0006 

0.0026 

-0.0628 

-0.0905 

-0.2381 

45 

0.0008 

0.0007 

0.0004 

-0.0407 

-0.0978 

-0.2354 

46 

0.0008 

0.0007 

0.0035 

0.1341 

-0.1020 

-0.2307 

47 

0.0007 

0.0010 

0.0075 

0.0162 

-0.1105 

-0.2277 

48 

0.0008 

0.0014 

0.0069 

-0.0321 

-0.1128 

-0.2212 

49 

0.0009 

0.0017 

0.0066 

0.0174 

-0.1143 

-0.2143 

50 

0.0010 

0.0020 

0.0053 

-0.0515 

-0.1149 

-0.2071 

MEANEI  through  MEANE4  represent  the  mean  of  estimation 
error  in  the  complete  state  vector  by  means  of  optimal 
filtering. 

MEANE1R  and  MEANE2R  represent  the  mean  of  estimation 
error  in  states  x,  and  x«  by  means  of  a  reduced  filter, 
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0.35 


1.5 
Time  in  seconds 

Fig.  3.3  -  Graphical  plot  of  the  variance  of  estima- 
tion error  in  xi  vs.  time.   Curves  a  and  b  show  the 
variance  of  estimation  error,  both  theoretical  and 
experimental,  for  the  optimal  estimation  process. 
Curves  c  and  d  show  the  corresponding  theoretical 
and  experimental  variances  for  the  reduced  case. 
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Fig.  3.4  -  Graphical  plot  of  the  variance  of  est. 
error  in  X2  vs.  time.   Curves  a  and  b  show  the  vari' 
ance  of  estimation  error,  both  theoretical  and  ex- 
perimental, for  the  optimal  estimation  process. 
Curves  c  and  d  show  the  corresponding  theoretical 
and  experimental  variances  for  the  reduced  case. 
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Fig.  3.5  -  Graphical  plot  of  the  means  of  estimation 
error  in  x^  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  Kalman  filter,  and 
curve  b  shows  that  for  the  reduced  filter. 
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Fig.  3.6  -  Graphical  plot  of  the  means  of  estimation 
error  in  X2  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  Kalman  filter,  and 
curve  b  shows  that  for  the  reduced  filter. 
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IV.  IMPLEMENTATION  OF  SUB-OPTIMAL  LINEAR 
FILTERS  FOR  DISCRETE  DYNAMIC  PROCESSES 
USING  MATRIX  PSEUDO-INVERSION 


In  the  previous  chapter  a  complete  discussion  of  the 
implementation  of  a  discrete,  reduced  model  Kalman  filter 
was  presented,  where  the  idea  of  partitioning  the  dynamics 
of  the  original  plant  was  principal.   In  this  chapter,  a 
different  approach  to  the  same  problem  is  presented, 
based  upon  the  work  of  Meditch  (Ref.3)  for  the  continuous 
case. 

Applying  the  same  idea  of  reduction  of  the  original 
process,  it  is  partitioned  for  this  case  into  two  or  more 
reduced  y.  vector  components  (i  =  1,2  ...  M) ,  with  each 
new  vector  containing  the  same  number  of  states.   This 
new  partitioning  is  accomplished  by  using  matrix  pseudo- 
inversion  (see  Appendix  A).   Once  the  original  system  has 
been  broken  down  into  M  reduced-order  sub-optimal  fil- 
ters, then  the  sub-optimal  estimates  for  the  original  pro- 
cess are  obtained  by  recombination  of  the  outputs  from  all 
the  sub-optimal  filters.   Section  A  presents  the  deriva- 
tion of  the  discrete  sub-optimal  filter  formulation. 
These  results  are  extended  to  Section  B,  where  expressions 


1 

The  order  (n)  of  the  original  system  is  equal  to  M 
times  the  order  (r)  of  the  reduced  plants,  that  is 
n  =  M  •  r. 
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for  the  sub-optimal  estimates,  theoretical  sub-optimal 
estimation  error  and  theoretical  sub-optimal  covariance  of 
estimation  error  are  derived.   Finally  in  Section  C  an  ex- 
ample is  solved  to  illustrate  the  technique. 

A.    DERIVATION  OF  THE  SUB-OPTIMAL  FILTER  EQUATIONS 

As  in  Chapter  III-A  it  is  assumed  that  the  system 
dynamics  are  represented  by  Eq .  (2.9),  without  consider- 
ing the  presence  of  the  deterministic  part  of  the  forcing 
vector,  that  is 

x(k+l)  =  $x(k)  +  Tw(k)  , 

with  measurement 

z(k)  =  Hx(k)  +  v(k)  , 

where  as  before,  w(k)  and  v(k)  are  random  gaussian  func- 
tions fully  described  by  Eqs.  (2.11)  through  (2.16). 
If,  as  discussed  in  Chapter  I,  it  is  desired  to 
estimate  only  r  states,  where  r^n ,  then  the  original  vec- 
tor x(k)  can  be  partitioned  or  transformed  into  two  or 
more  vectors  by  using  a  set  of  constant  matrices 

Ai,    i  =  1,2  . ..  M 

with  A.  having  dimensions  (r  x  n).  Notice  that  each  A. 
must  have  r  rows  and  n  columns.  Then  it  is  possible  to 
define  a  new  set  of  reduced  vectors 
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w 


Zi(k)  =  A.x(k),    i  =  1,2  ...  M  (4.1) 

ith  y_.(k)  representing  a  (r  x  1)  vector. 

Now  the  concept  of  the  Generalized  Inverse   of  a 

matrix  (pseudo-inversion)  is  applied  to  the  implementa- 

+  2 
tion  of  the  sub-optimal  filters.   A  matrix  A    is  con- 
sidered to  be  the  pseudo-inverse  of  a  rectangular  matrix 
A  if 


AA+A  =  A  ,  (4.2) 


where  in  general  A  is  a  (s  x  t)  matrix  and  A   is  a  (t  x  s) 

matrix. 

The  selection  of  the  set  of  constant  matrices  A.  will 

i 

have  as  a  necessary  requirement  for  the  solution  of  the 
filter  implementation  and  the  proper  recovery  of  the  sub- 
optimal  estimates  from  the  combination  of  the  filter  out- 
puts, that 


M 

E 

i=l 


AtA.  =  I  ,  (4.3) 

i  l      ' 


1 

See  Deutsch  (Ref.10)  and  Penrose  (Ref.8)  for  a  complete 
discussion  about  matrix  pseudo-inversion. 


A   stands  for  the  pseudo-inverse  of  A. 
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where   I   is  the  identity  matrix.   This  requirement  will 
be  clarified  after  the  derivation  which  follows. 

If  Eq .  (4.1)  is  premult iplied  by  A. ,   i  =  1,2  ...  M 
and  summation  over  all  the  ifs  is  performed,  then 


M  M 

V  A+y.(k)  =  T     A+A.x(k) 
i=l  i=l 


from  which  it  is  clearly  seen,  applying  the  requirement 
established  by  Eq .  (4.3),  that 


M 

x(k)  =     A+y_.(k)  .  (4.4) 

i=l 


This  last  equation  provides  the  key  for  recovery  of 
the  original  system  state  vector. 

Now  if  the  original  discrete  system  equation  is  pre- 
multiplied  by  A. ,   i  =  1,2  . . .  M, 

Aix(k+1)  =  A^xCk)  +  ATw(k)  ,  (4.5) 

and  if  Eqs.  (4.1)  and  (4.4)  are  substituted  in  Eq .  (4.5), 
the  latter  becomes 

M 
y^Ck+l)  =  Ai$  ^T  A^y_j(k)+Airw(k),   i  =  1,2  ...  M 
j=l 

(4.6) 
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Expression  (4.6)  represents  a  set  of   M   reduced,   dis- 
crete system  equations,  each  of  order  r. 

The  measurement  process  can  be  written  assuming  that 
all  the  measurements  that  were  available  for  the  complete 
filter  are  available  too  for  the  reduced  filters. 

z(k)  =  Hx(k)  +  v(k)  , 

becomes  after  replacement  of  x(k)  by  Eq .  (4.4), 


z(k)  = 


M 

h  y  Atx  (k)  +  v(k) 

3=1 


(.47) 


Fig.  (4.1)  will  help  to  better  visualize  the  partitioned 
system. 

In  order  to  facilitate  the  following  derivations,  the 
optimal  Kalman  filter  estimation  Eq .  (3.8)  may  be  rewritten 
in  a  slightly  different  form,  that  is 


x(k/k)  =  $x(k-l/k-l)  +  G(k)[z(k)  -  z(k/k-l)J 


(4.8) 


where  £(k/k-l)  represents  the  expected  value  of  ^(k)  given 
all  measurements  up  to  time  (k-1),  or 


|(k/k-l)  =  E|z(k)|     given   z(k-l)  . 
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From  Eq.  (4.7),  however, 


M 


;(k)  =  H  Y,    Aj  Xj<k)  +  vCk) 
j-l 


(4.7) 


Therefore , 


z(k/k-l)  = 


M 

H  Y  At  E'y,(k> 

' — '    J    L  J 

j-l 


given  z-(k-l) 


or 


M 


z(k/k-l)  = 


H  [  A+  ^(k/k-1) 
j-l 


(4.9) 


Considering  further  that 


x(k/k-l)  =  Jx(k-l/k-l)  , 


(4.10) 


where  x(k/k-l)  represents  the  expected  value  of  x(k)  given 
all  measurements  up  to  time  (k-1),  and  premult iplying  Eq . 
(4.10)  by  At,  1,2  ...  M, 


Aix(k/k-l)  =  Ai^x(k-l/k-l)  . 


(4.11) 


By  extension  of  Eq.  (4.1) 


A.£(k/k-l)  =  2i(k/k-l) 


49 


and  similarly  by  proper  use  of  Eq .  (4.4) 


x(k-l/k-l)  = 


M 

£  A+  £.(k-l/k-l)  , 

j-l 


so  Eq .  (4.11)  becomes 


M 


£i(k/k-l)  =  kti    £  At ^(k-iyk-i) 

j-l 


(4.12) 


Basing  the  implementation  of  the  reduced  filters  upon 
the  dynamics  of  the  set  of  systems  described  by  Eq.  (4.6), 
a  set  of  M  reduced-filter  estimation  equations  may  be 
stated  by  induction  from  Eq.  (4.8),  as 


y_.(k/k)  =  ^(k/k-l)  +  G.(k)[z(k)  -  |(k/k-l) 


Replacing  for  ^(k/k-1)  and  z(k/k-l)  Eqs.  (4.12)  and  (4.9), 
gives 

M  M 

2i(k/k)=Ai$  Y,     A+£j(k-l/k-l)+G1(k)[z(k)-H  V  A^k/k-l)] 


(4.13) 


with  i  =  1,2  ...  M. 


th 


This  expression  represents  a  set  of  M  r   -order  coupled 
reduced  filters,  whose  outputs  £.  (k/k)  will  be  recombined 
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in  a  proper  way  in  order  to  get  the  sub-optimal  estimates 
x(k/k),  where  x(k/k)  is  an  (n  x  1)  vector  (see  Section 
IV-B). 

Referring  to  the  optimal  gain  and  theoretical  vari- 
ance expressions  for  the  Kalman  filter  (Eqs.  3.7a,b,c), 
the  equivalent  M  subsets  of  sub-optimal  equations  for  each 
one  of  the  reduced  filters  may  be  expressed,  by  definition, 
as 

G.(k)  =  Pi(k/k-l)H^[HiPi(k/k-l)H^+R]    ,       (4.14a) 


where 


+ 


H._  =  H  AT    (see  Eq .  4.7) 


corresponds  to  the  i    subsystem  contribution  to  the 
measurement  vector; 


Pi(k/k)  =[l  -  Gi(k)Hjpi(k/k-l)  , 


(4.14b) 


and 


Pi(k+l/k)  =  ^P^k/k)^  +  Q 


(4.14c) 


where 


(jh  =  A^At    (see  Eq .  4.6) 
Qi  =  A.rE[w(k)wT(k)]rV 


T 

a.qa! 

lx  1 
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and 


r.  =  a  .  r  . 
i   i 


It  is  important  to  notice  that  although  the  M  estima- 
tion equations  and  the  measurement  equation  are  coupled 
for  all  the  M  reduced  filters,  the  gain  and  covariance 
equations  for  the  i    filter  are  truly  uncoupled  from  the 
other  j  =  1,2  ...  M  filters,  where  j  ^   i. 

Eqs.  (4.14)  complete  the  implementation  of  the  M  sub- 
optimal  filters.  The  filtering  process  will  be  visualized 
better  by  means  of  Fig.  (4.2). 

The  initialization  values  for  the  sub-optimal  filters 
are  derived  from  the  statistics  used  to  initialize  the  com- 
plete optimal  filter  in  the  following  way.   Given  the 
values  for  P(0/-1)  and  knowing  that 


P.(0/-1)  =  E  ^(0)^(0) 


T, 


and 


Xt(0)  =  Aix(0)  , 


then 


P.(0/-1)  =  EA,x(0)xT(0)aT1 

1  L    1 1 


A.E[x(0)xT(0)]A^ 


A.P(0/-1)A^  . 


(4.15) 
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By  the  same  reasoning, 

£.(0/-l)  =  A.£(0/-1)  (4.16) 

will  be  used  as  the  initialization  values  for  the  set  of  M 
reduced  estimation  Eqs.  (4.13). 

B.    SUB-OPTIMAL  ESTIMATES,  ESTIMATION  ERROR  AND  SUB- 
OPTIMAL  VARIANCE  OF  ESTIMATION  ERROR 

In  the  preceding  section,  M  reduced  filters  were  im- 
plemented after  breaking  down  the  original  process  into  M 
reduced  (r  x  l)  vectors.  In  order  to  obtain  the  sub-optim- 
al estimates  of  the  original  process,  the  outputs  y. (k/k) 
from  the  reduced  filters  have  to  be  recombined.  This  is 
easily  accomplished  by  proper  application  of  Eq .  (4.4)  to 
the  filter's  estimates,  whereupon 

M 
|(k/k)  =  £  A+  £.(k/k)  (4.17) 

j-i 

will  represent  the  sub-optimal  estimates  of  the  original 
process,  with  x(k/k)  representing  an  (n  x  1)  vector. 

The  remainder  of  this  section  is  devoted  to  the  de- 
velopment of  a  matrix  difference  equation  for  the  itera- 
tive computation  of  the  sub-optimal  covariance  of  estima- 
tion error. 

If  the  error  of  sub-optimal  estimation  is  defined  as 

e(k)  =  x(k)  -  x(k/k)  (4.18a) 
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then  using  Eq .  (4.17)  it  becomes 

M 

e(k)  =  x(k)  -  V  At  £.(k/k)  , 
J-l     J 


(4.18b) 


or 


M 


e(k+l)  =  x(k+l)  -  V    At  £.(k+l/k+l)  . 

— ■   j   j 

j-l 


(4.18c) 


Substitution  of  Eqs.  (4.13)  and  (4.7)  into  Eq .  (4.18c) 
gives 

M     f     M 

e(k+l)  =  x(k+l)  "  J]    At  Ai$  Z  At  2*(kA) 

— '   J  L     — '   J   J 

j-l         j-l 

M  M  1 

+Gi(k+1)[H     £    A+  Zj(k+l)+v(k+l)-H$     £    Aj2j(kA)]J 
j-l  j-l 

M 
Replacing       )      A.y_.(k+1)    with  x(k+l),    and    substituting 

j-l 

x(k+l)    =   $x(k)   +   Tw(k)    from   Eq .    (3.1),    then 


M 


e(k+l)    =   $x(k)    +   rw(k)    - 


At 


j-i 


M 
A.$     V     A+£.(k/k) 

j-l 


M 


+Gi(k+l)[H$x(k)    +   Hrw(k)    +   v(k+l)    - 


H$    £    Aj2j(kA) 
j=l 
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Expansion  of  this  expression  for  the  sub-optimal  error 
leads  to 

M       :m 

e(k+l)  =  $x(k)  +  rw(k)  -  V     A"Ja.$  )   A+£.(k/k) 
-         —  ' —   J  J   L — '   J  J 

j-l        j=l 


M 


5  A"!"G.(k+l)H$x(k) 
j-l 


M 

Y  A"t"G.(k+l)Hrw(k) 
j-l 


M  KM 

WtG.(k+l)v(k+l)+  Y,A'tG.(k+l)H$  y  A+z-k/k). 

j-l  j-l  j-l 


Collecting  terms, 


M  M 

e(k+l)  =  [\   -  £  A+G.(k+l)H$][x(k)  -  ^A+^Ck/k)] 
j=l  j=l 


M 


+  !I-   UtG.(k+l)H  ]  Tw(k)- 
j=l 


M 

VA^G^+lMk+l) 
j=l 


Substitution  of  Eq .  (4.18b)  for  e_(k)  gives 


M 


M 


e(k+l)  = 


+, 


I-   ,  ATG.(k+l)H'^e(k)+ 
j=l  J  J 
M 

Y  A"!"G.(k+l)v(k+l). 
j=l  J  J 


I  - 


V  AtG.(k+l)H]rw(k) 
j=l  J  J 

(4.19) 
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Eq .  (4.19)  represents  a  theoretical  form  for  the  error 
of  sub-optimal  estimation  which  can  be  useful  in  establish- 
ing the  steady-state  estimation  error  for  the  sub-optimal 
filter. 

If  the  sub-optimal  theoretical  variance  of  estimation 
error  is  defined  as 


S(k)  =  E!e(k)er(k)J 


(4.20) 


where  ^(k)  is  an  (n  x  n)matrix,  then  substitution  of  Eq . 
(4.19)  for  eXk)  gives 


S(k+1)  =  E\ 


M 


I  -  ]T  AtG.(k+l)H](j)e(k)  + 
j-l 


M 
I-  Y  AtG.(k+l)H|rw(k) 
j  =  l 


M 
-  )^A"t"G.(k+l)v(k+l) 
j-l 


same  terms 


T 


Carrying  out  the  indicated  matrix  operations  term  by  term 
yields 


M 


M 


5(k+l)  =  E<[l-  y  AtG.(k+l)H]$e"(k)e"T(k)$T("l-  V  A"I"G . (k+1) H 


.IT 


M 


M 


+ 


I-  y  AtG.(k+l)Hlrw(k)wT(k)rT|"l-  )  AtG.(kfl)H 

J-l  j-l 


T 


M 


M  1 

+  ^A+Gj(k+l)v(k+l)yT(k+l)  [G!(k+1)A+T+  products 
j-l  j-l  J 
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Given  the  statistics  of  the  problem,  all  the  cross  products 
will  yield  zero  expectation.   S(k+1)  then  becomes 


M 


M 


5(k+l)  =  [i  -  ]T  AJG^k+DHJpfeCk^Ck)]  $T  [i  -  Y  AtG^k+DH 


-.T 


j  =  l 


j-l 


M  M 

+  [l  -    ^  AtG.(k+l)H]rE[w(k)wT(k)]rT[l-  VAtc(k+l)H 
j  =1    J    J  j = 1    J    J 


,T 


M 


M 


+  )    AtG.(k+l)E[v(k+l)vT(k+l)]    V  GT(k+l)A"tT    , 


j=l 


j-l 


and    finally 


M 

S(k+1)    =    !I    -  V  Atc(k+1)H   (j>S(k)(|)T 
j  =  l 


M 


I    - 


+. 


)    Ap.(k+1)H 

j-i 


T 


M 


+ 


+  I    -    )    A"G.(k+l)H!Q  |I    - 


j=l 


J    J 


m 


M 

A"!"G.(k+l)H 
J    J 

j-l 


T 


M 


M 


+  AtG.(k+l)R  GT(k+l)A"l"T  . 

—     J    J  J  J 

j-l  j-l 


(4.21) 


Eq .    (4.21)    represents    an    iterative   expression    for   the 
computation    of    the    theoretical    sub-optimal   covariance    of 


58 


estimation  error.   To  initialize  the  computation  of  SXk+l), 
the  initial  value  of  3"(0)  should  be  chosen  as 


3(0)  =  E  e(0)e~T(0) 


M  2 

=  E  x(o)  -    a1"£.(o/-i)  ; 


2 
=  E  Lx(0)  -  x(0/-l) 


=  P(0/-1)  . 

The  discrete  matrix  expression  for  the  theoretical 
sub-optimal  covariance  of  estimation  error  3(k),  may  be 
used  as  a  reference  to  evaluate  the  performance  of  the  re- 
duced set  of  filters,  i.e.,  how  close  to  the  "true"  values 

are  the  sub-optimal  estimates.   The  solution  for  S"(k) 
should  be  compared  with  the  experimental  variance  of  sub- 
optimal  estimates. 

C.    EXAMPLE  2 

The  model  used  to  investigate  and  to  illustrate  the 
technique  presented  in  this  chapter,  was  exactly  the  same 
fourth-order  plant  used  in  Example  1,  subjected  to  the 
same  deterministic  and  random  forcing  functions  and  with 
precisely  the  same  measurements.   The  continuous  system 
was  discretized  using  the  same  sampling  time  T  =  0.05  sec, 
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and  therefore  the  same  $  and  T    matrices  were  obtained, 
arriving  at  the  discretized  system  described  by 

x(k+l)  =  $x(k)  +  r   u(k)  +  w(k)   ,  (2.7) 

with  measurement 

z(k)  =  Hx(k)  +  v(k)  . 

The  goal  of  this  example,  as  in  Example  1,  was  to 
estimate  only  two  states  of  the  fourth-order  system 
(x-j  and  X2)  by  means  of  a  reduced-order  filter.   Apply- 
ing the  technique  discussed  in  Section  A  of  this  chapter, 
the  original  model  (4   order)  could  be  partitioned  into 
two  second-order  models,  but  in  order  to  make  the  compu- 
tational procedures  much  simpler,  the  original  model  was 
partitioned  into  four  scalar  models.   This  partitioning 
approach  facilitated  the  selection  of  the  A.  matrices  and 
the  computation  of  their  respective  pseudo-inverses,  A. . 

The  solution  to  the  partitioning  problem  was  achiev- 
ed by  trial  and  error.   Four  new  vectors  y_.  ,  i  =  1,2,3,4 
were  defined  by 

Yi(k)  =  Aix(k)  ,   i  =  1,2,3,4. 


The  four  A.  (1  x  4)  matrices  selected  were 
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A,  =   1   0   0   0 

1         L 


A2  =  :  0   1   0   0 


A3  =   0   0   1   0 


A,  =      0   0   0   1 

4    < 


It  is  important  to  remember  here  that  all  the  A.  matrices 

should  have  the  same  dimensions,  and  that  the  choice  of  the 

A. 's  in  this  case  renders  the  y. 's  scalar  quantities. 
i  J  l  ^ 

The  selection  of  the  A.  matrices  was  made  by  trial  and 
error  (in  this  case  obvious  and  simple)  in  order  to  get  the 
values  of  proper  A. fs   that  would  fulfill  the  requirement 
established  by  Eq .  (4.3),  in  this  particular  case 


4 

V 

<L_ 
1=1 


V  A+A.  .  I. 

^  ,   Li 


The  pseudo-inverse  for  each  A.  matrix  was  found  using 
Eq.  (4.2) 


A.  At  A,  =  A,  ,  (4.2) 


1 

Werther  (Ref.9)  presents  a  computation  method  that 
yields  to  matrix  inversion  or  pseudo-inversion. 
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from  which  the  four  A.  (4  x  1)  matrices  were  found  to  be 


4- 


1 
0 
0 
0 


4  = 


o 
l 

0 
0 


At  = 


0 
0 

1 


A   = 


rol 


The  partitioning  of  the  original  state  vector  is  then 
summarized  by 


*1 

^4 


=  : 1   0   0   0  x(k)  = 


xx(k) 


[0   1   0  Ojx(k)  =  x2(k) 
JO   0   1   Ojx(k) 

[o   0  0   lix(k) 


xQ(k) 


x4(k) 


with  dynamics  as  specified  by  Eq .  (4.6) 

This  partitioning  can  be  interpreted  as  an  apparent 
canonization  of  the  fourth-order  model  into  four  first- 
order  models,  represented  each  one  by  one  state  of  the 
original  system. 

Having  selected  the  A.  and  A.  matrices,  four 

&  i      i         ' 

coupled  scalar  filters  were  implemented  in  parallel  using 


It  can  be  easily  verified,  that  the  selection  of  the  A. 
and  A.  matrices  agrees  with  Eqs.  (4.2)  and  (4.3). 
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Eqs.  (4.13)  and  (4.14),  but  due  to  the  presence  of  a  deter- 
ministic forcing  function,  Eq .  (4.13)  was  changed  slightly 
as  follows, 


Yi(k/k) 


£  +  r 

=  A.$  >_  A^yi(k-l/k-l)  +  nu(k-l)+Gi(k)|z(k) 


j-l 


-H   Aj[yj(k/k-l)  +  r.u(k-l)]   1  =  1,2,3,4 
j=l 

(4.22) 

to  show  the  effect  of  u(k)  on  the  prediction  step. 

The  following  dynamics  were  used  for  the  computation 
of  the  filter  gains  and  theoretical  covariances. 


+ 


$.  =  a.  (5  a? 
*i   1*  i 


i  =  1,2,3,4 


1   0   0   0 


h  = 


11 


22 


41 


'33 


14 


i  44 


1 

0 
0 
0 


fl     '11  ' 
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where  './••-,  and  /.^  are  scalars.   Similarly, 


i    r 


[O   1   O   Oj 


•/" 


/ 
I  2 


22 


i 


33 


L  '41 


14 

roi 

1 

i 

i  o 

44  _ 

0 

1.2     '  A.22.  ,' 

'3  =  '33  ' 


and 


/4  ~  /44  • 


Also, 


r.  = 


A.r 

i 


i  =  1,2,3,4 


where 


7,  = 

"  1 


i  o  o  o]|"7nl 

ki! 
j 

7! 
31 

'41 


-ru 
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and  similarly, 


72  ~  72l 

'h  =  hi 


and 


7  -T,i 

"4    ''41 


Note  that  the  symbols  yU  ,  'f\  >  •••   -i  >  '9  •  •  •  etc«> 
correspond  to  scalar  dynamics  of  the  reduced  processes, 

whereas  the  symbols  /it>  /o?'  ' "'  7ll»  7?1'  *'*  etc*> 
correspond  to  scalar  elements  of  the  dynamics  of  the  full 
system. 

The  Q  matrix  was  partitioned  as  follows: 


Q. 

l 


A.QAt 
i^  l 


i  =  1,2,3,4 


then 


I' 


10  0  0 


11 


'41 


22 


33 


14 


44 


1 
0 
0 

0 


H 


n 
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and  similarly 


q2    q22  ' 
q3    q33    ' 


and 


q4  "  q44  ' 

where  q-.  ,  qo  •••  q-i  -i  >  q??  •••  etc«>  represent  scalar  quan- 
t it  ies . 

The  measurement  process  was  specified  by 


+ 


z(k)  =  H      AjYj (k)  +  v(k) 
3=1 


with  the  measurement  matrix  H  =   1   0   0   0  j ,  and 
v(k)  :  N(0.0,  0.25). 

For  gains  and  variance  calculations  for  each  filter, 
the  measurement  matrix  H(l  x  4)  was  partitioned  as  fol- 
lows : 


H.  =  HAt 
l      l 


i  =  1,2,3,4 


from  which 


hl  " 


1   0   0   0  I    1 


0 

L  o  J 
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or 


h1  =  1, 
h2  =  0, 

h3  -  0, 


and 


\   =   0, 


with  h-,  ,  h.2  ...  being  scalars. 

The  initialization  values  for  the  reduced  filters 
were  computed  using  Eqs.  (4.15)  and  (4.16),  from  which 
the  initial  guesses  for  the  covariance  of  estimation 
error  P.(0/-1)  and  expected  value  of  the  state  vector 
y.(0/-l)  were  obtained  as  follows: 


P.(0/-1)  =  A.P(0/-1)A+      i  =  1,2,3,4 


which  gives 

[l  o  o  o]rPll(o/-i)  |[  1 

p22(0/-1) 

p33(o/-i) 


P^O/-!)  = 


L 


P44(0/-1) 


0 

i 


0 
0 
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or 


Pl(0/-1)  =  pn(0/-l)  =  0.5  , 
p2(0/-l)  =  p22(0/-l)  =  2.0  , 
p3(0/-l)  =  p33(0/-l)  =  2.0  , 


and 


p4(0/-l)  =  p44(0/-l)  =  4.0  , 

where  p-,(0/-l)  ...,  p-.-i(0/-l)  ...,  are  scalars,  with 
p-,(0/-l),  P2(0/-1),  ...  etc.,  representing  the  initial 
values  of  covariance  for  the  reduced  filters. 

For  the  initial  values  of  the  estimates  it  is  easily 
seen  that 

y^OAl)  =  ^(0/-l)  =  0.0  , 
y2(0/-l)  =  x2(0/-l)  =  0.0  , 
y3(0/-l)  =  x3(0/-l)  =  0.0  , 

and 

y4(0/-l)  =  x4(0/-l)  =  0.0 

will  provide  the  filters  with  unbiased  initial  values. 
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Then  four  scalar,  discrete,  sub- optimal  filters,  in 
this  example  coupled   only  through  the  measurement  process, 
were  implemented  and  run  in  parallel  in  accordance  with 
Eqs.  (4.22),  (4.14),  and  (4.7).   At  the  end  of  each  sam- 
pling interval  the  outputs  from  the  four  filters  were  com- 
bined by  the  expression 

4 
x(k/k)  =  V     A+y\(k/k)  , 

j-i 

< 

in  order  to  recover  the  sub-optimal  estimates  of  the 
original  vector  x(k),  from  which  only  the  estimation  of 
Xi  and  Xj   was  of  particular  interest. 

At  the  same  time,  the  theoretical  sub-optimal  covari- 
ance  of  estimation  error  ^"(k)  ,  was  computed  using  recur- 
sive matrix  Eq.  (4.21),  in  order  to  evaluate  the  selection 
of  the  set  of  constant  matrices  A. .   This  is  possible  to 
see  when  S~(k)  is  compared  with  the  experimental  reduced 

covariance  of  estimation  error. 

2 
Again  a  Monte  Carlo  simulation   was  performed  with 

200  replications  of  the  trajectory,  each  trajectory  made 


1 

This  is  not  generally  so.   The  uncoupling  here  is  a  con- 
sequence of  the  choice  of  the  Aj_  matrices.   The  y]_    state 
equations  need  not  be  uncoupled,  but  the  associated  gain 
and  variance  equations  are,  however,  always  uncoupled. 

2 

As  in  Example  1,  a  computer  System  360/67  was  used  for 
the  simulation  and  solution  of  this  problem.   See  Appendix 
C  for  corresponding  computer  program. 
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of  50  samples,  with  sampling  intervals  of  0.05  seconds.  As 
a  result  of  this  simulation,  ensemble  values  for  the  mean 
and  experimental  covariance  of  estimation  error,  as  a  func- 
tion of  time,  were  obtained  in  x-,  and  X2  for  the  optimal 
and  sub-optimal  filters  as  a  means  of  comparing  their  per- 
formance. 

Numerical  results  of  this  simulation  are  tabulated  in 
Tables  (IV-2)  through  (IV-6)  and  graphical  representation 
of  the  results  is  shown  in  Figs.  (4.3)  through  (4.6). 
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Table 

IV-1 

PHI  MATRIX 

0.9999 

0.0499 

0.0012 

0.0000 

-0.0066 

0.9946 

0.0480 

0.0009 

-0.3644 

-0.2981 

0.8871 

0.0307 

-12.2713 

-10.1815 

-3.9181 

0.3043 

GAMMA  MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 

Q  MATRIX 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0008 

0.0000 

0.0000 

0.0013 

0.0447 

0.0000 

0.0008 

0.0447 

1.5059 

VARIANCE  OF  FORCING  NOISE  =0.01 
VARIANCE  OF  MEASUREMENT  NOISE  =0.25 

INITIALIZATION  -  P(0/-1) 


0.5000 

0.0 

0.0 

0.0 

0.0 

2.0000 

0.0 

0.0 

0.0 

0.0 

2.0000 

0.0 

0.0 

0.0 

0.0 

4.0000 

INITIALIZATION  -  X(0/-l) 

0.0 
0.0 
0.0 
0.0 
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Table   IV-2 
GAINS    FOR    FULL    AND   RFOUCED    FILTERS 


G(l»l) 


G(2,l) 


GOtll 


GUtl)         Gllhll         G2(l,l) 


1 

0.6667 

2 

0.4071 

3 

0.3097 

4 

0.2677 

5 

0.2477 

6 

0.2349 

7 

0.2227 

8 

0.2087 

9 

0.1928 

10 

0.1758 

11 

0.1584 

12 

0.1412 

13 

0.1247 

14 

0.1093 

15 

0.0950 

16 

0.0819 

17 

0.0701 

18 

0.0596 

19 

0.0502 

20 

0.0420 

21 

0.0348 

22 

0.0286 

23 

0.0233 

24 

0.0188 

25 

0.0151 

26 

0.0120 

27 

0.0095 

28 

0.0075 

29 

0.0059 

30 

0.0047 

31 

0.0037 

32 

0.0031 

33 

0.0026 

34 

0.0022 

35 

0.0020 

36 

0.0018 

37 

0.0018 

38 

0.0017 

39 

0.0017 

40 

0.0017 

41 

0.0017 

42 

0.0017 

43 

0.0017 

44 

0.0017 

45 

0.0017 

46 

0.0017 

47 

0.0017 

48 

0.0017 

49 

0.0017 

50 

0.0017 

0.0 

0.2333 

0.4075 

0.5088 

0.5314 

0.4901 

0.41Q8 

0.3175 

0.2268 

0.1474 

0.0825 

0.0322 

-0.0051 

-0.0315 

-0.0490 

-0.0596 

•0.0647 

-0.0660 

-0.0643 

-0.0606 

-0.0556 

-0.0499 

-0.0438 

-0.0378 

-0.0320 

-0.0266 

-0.0217 

-0.0174 

-0.0136 

-0.0105 

-0.0078 

■0.0057 

•0.0041 

-0.0028 

-0.0018 

-0.0011 

•0.0006 

•0.0003 

-0.0001 

0.0000 

0.0001 

0.0001 

0.0001 

0.0001 

0.0000 

0.0000 

■0.0000 

•0.0000 

•0.0001 

•0.0001 


0.0 

-0.2095 

■0.6537 

•1.1719 

•1.6080 

•1.8640 

-1.9266 

•1.8378 

-1.6557 

-1.4305 

-1.1970 

•0.9760 

■0.7780 

•0.6072 

•0.4637 

-0.3457 

■0.2505 

-0.1752 

•0.1167 

■0.0723 

-0.0395 

■0.0161 

-0.0001 

0.0100 

0.0157 

0.0181 

0.0181 

0.0166 

0.0141 

0.0111 

0.0080 

0.0050 

0.0023 

•0.0001 

■0.0022 

•0.0038 

•0.0050 

•0.0060 

•0.0066 

•0.0070 

•0.0073 

•0.0074 

•0.0074 

•0.0073 

-0.0072 

•0.0071 

•0.0070 

■0.0068 

■0.0067 

■0.0066 


0.0 

-7.2845 

■9.4633 

■9.0573 

•7.0288 

■4.4015 

•1.9563 

-0.0721 

1.1851 

1.9077 

2.2348 

2.2961 

2.1931 

1.9984 

1.7610 

1.5120 

1.2705 

1.0474 

0.8480 

0.6746 

0.5269 

0.4036 

0.3026 

0.2215 

0.1576 

0.1084 

0.0714 

0.0444 

0.0254 

0.0127 

0.0046 

■0.0000 

•0.0022 

•0.0028 

•0.0024 

•0.0015 

■0.0004 

0.0007 

0.0017 

0.0024 

0.0030 

0.0033 

0.0034 

0.0033 

0.0031 

0.0028 

0.0025 

0.0022 

0.0018 

0.0015 


0.6667 
0.4000 
0.2857 
0.2222 
0.1817 
0.1538 
0.1333 
0.1176 
0.1052 
0.0952 
0.0869 
0.0799 
0.0740 
0.0689 
0.0644 
0.0605 
0.0571 
0.0540 
0.0512 
0.0487 
0.0464 
0.0444 
0.0425 
0.0407 
0.0391 
0.0377 
0.0363 
0.0350 
0.0338 
0.0327 
0.0317 
0.0307 
0. 0298 
0.0289 
0.0281 
0.0273 
0.0266 
0.0259 
0.0252 
0.0246 
0.0240 
0.0234 
0.0229 
0.0224 
0.0219 
0.0214 
0.0210 
0.0205 
0.0201 
0.0197 


0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 


G  (4x1)    represents  the  matrix  of  gains  for  the  optimal 
filter. 

Gl  and  G2  represent  the  gains  used  in  the  sub-optimal 
filters  1  and  2. 
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Table    IV-3 


THEORETICAL    VARIANCE    OF    ESTIMATION    ERROR 
FULL     PLANT 


K 

PK( 1, i) 

PK(2,2I 

PK<3,3) 

PK<4,4) 

1 

0.1667 

2.0000 

2.0000 

4.0000 

2 

0,1018 

1.9603 

1.7605 

242.6279 

3 

0.0774 

1.8002 

2.3189 

346.3362 

4 

0.0669 

1.5182 

3.9792 

261.1479 

5 

0.0619 

1.1771 

5.7036 

142.2309 

6 

0.0587 

0.8478 

6.6569 

60.6997 

7 

0.0557 

0.5745 

6.6825 

21.3200 

8 

0.0522 

0.3709 

6.0419 

8.2223 

9 

0.0482 

0.2307 

5.0782 

6.9855 

10 

0.0439 

0.1400 

4.0514 

0.3412 

11 

0.0396 

0.0847 

3.1097 

11.6978 

12 

0.0353 

0.0529 

2.3160 

12.9672 

13 

0.0312 

0.0361 

1.6824 

13.1242 

14 

0.0273 

0.0282 

1.1960 

12.477R 

15 

0.0237 

0.0255 

0.8337 

11.3650 

16 

0.0205 

0.0254 

0.5710 

10.0516 

17 

0.0175 

0.0264 

0.3855 

8.7171 

18 

0.0149 

0.0275 

0.2581 

7.4696 

19 

0.0126 

0.0283 

0.1736 

6.3652 

20 

0.0105 

0.0286 

0.1199 

5.4260 

21 

0.0087 

0.0282 

0.0877 

4.6526 

22 

0.0072 

0.0273 

0.0701 

4.0336 

23 

0.0058 

0.0259 

0.0619 

3.5512 

24 

0.0047 

0.0241 

0.0595 

3.1854 

25 

0.0038 

0.0221 

0.0601 

2.9159 

26 

0.0030 

0.0199 

0.0619 

2.7237 

27 

0.0024 

0.0177 

0.0639 

2.5918 

28 

0.0019 

0.0155 

0.0654 

2.5057 

29 

0.0015 

0.0133 

0.0659 

2.4529 

30 

0.0012 

0.0114 

0.0654 

2.4236 

31 

0.0009 

0.0096 

0.0640 

2.4099 

32 

0.0008 

0.0081 

0.0618 

2.4059 

33 

0.0006 

0.0067 

C.0590 

2.4071 

34 

0.  0006 

0.0056 

0.0557 

2.4105 

35 

0.0005 

0.0046 

0.0523 

2.4141 

36 

0.0005 

0.0039 

0.0489 

2.4167 

37 

0.0004 

0.0033 

0.0455 

2.4177 

38 

0.0004 

0.0028 

0.0424 

2.4170 

39 

0.0004 

0.0024 

0.0396 

2.4146 

40 

0.0004 

0.0022 

0.0371 

2.4109 

41 

0.0004 

0.0020 

0.0349 

2.4062 

42 

0.0004 

0.0018 

0.0330 

2.4008 

43 

0.0004 

0.0017 

0.0315 

2.3951 

44 

0.0004 

0.0017 

0.0303 

2.3894 

45 

0.0004 

0.0017 

0.0293 

2.3839 

46 

0.0004 

0.0016 

0.0285 

2.3787 

47 

0.0004 

0.0016 

0.0279 

2.3741 

48 

0.0004 

0.0016 

0.0275 

2.3700 

49 

0.0004 

0.0016 

0.0272 

2.3664 

50 

0.0004 

0.0016 

0.0270 

2.3635 

PK  represents   the  diagonal   elements   of   the   theoretical 
covariance  matrix   for   the   optimal    filter. 
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Table  IV-4 


EXPERIMENTAL  VARIANCE  OF  ESTIMATION  ERROR 
FULL  PLANT 


K 

VK( 1,1) 

VM2,2) 

VM3,3> 

VM4,4) 

1 

0.1317 

1.9066 

1.9733 

4.0880 

2 

0.0910 

1.8538 

1.7402 

225.0217 

3 

0.0739 

1.7340 

2.2403 

326.4661 

4 

0.0623 

1.4891 
1.1562 

3.7969 

257.7678 

5 

0.0645 

5.6106 

138.7382 

6 

0.0518 

0.8086 

5.9918 

60.5067 

7 

0.0510 

0.5760 

6.5021 

23.3325 

8 

0.0485 

0.3807 

6.0129 

8.6785 

9 

0.0465 

0.2391 

5.1270 

7.2472 

10 

0.0459 

0.1500 

4.3535 

9.5983 

11 

0.0424 

0.0927 

3.5224 

12.3814 

12 

0.0394 

0.0557 

2.7024 

14.0400 

13 

0.0367 

0.0360 

1.9633 

13.9882 

14 

0.0295 

0.0276 

1.3043 

13.4596 

15 

0.0264 

0.0240 

0.9568 

11.7134 

16 

0.0225 
0.0187 

0.0239 

0.6456 

12.1689 

17 

0.0256 

0.4017 

9.6542 

IS 

0.0152 
0.0130 

0.0259 
0.0271 

0.2554 

7.2280 
5.9514 

0.1698 

20 

0.0111 

0.0269 

0.1194 

5.5124 

\\ 

0.0093 

0.0275 

0.0835 

5.0569 

0.0073 

0.0253 

0.0671 

4.1066 

23 

0.0060 

0.0241 

0.0626 

3.5153 

24 

0.0047 

0.0218 

0.0616 

2.8292 

25 

0.0039 

0.0203 

0.0586 

3.1777 

26 

0.0031 

0.0186 

0.0564 

2.4777 

27 

0.0025 

0.0166 

0.0546 

2.7436 

28 

0.0020 

0.0150 
0.0131 

0.0546 

2.7655 

29 

0.0016 

0.0550 

2.6672 

30 

0.0012 

0.0115 

0.0581 

2.0752 

31 

0.0010 

0.0099 

0.0586 

2.4300 

32 

0.0008 

0.0084 

0.0572 

2.7947 

33 

0.0007 

0.0072 

0.0546 

2.5196 

34 

0.0006 

0.0060 

0.0495 

2.1787 

35 

0.0005 

0.0051 

0.0437 

2.5699 

36 

0.  0004 

0.0043 

0.0422 

2.2536 

37 

0.0004 

0.0035 

0.0438 

2.5338 

38 

0. 0004 

0.0029 

0.0449 

2.3031 

39 

0.0004 

0.0025 

0.0427 

2.0479 

40 

0.0004 

0.0021 

0.0373 

2.4109 

41 

0.0004 

0.0019 

0.0308 

2.2848 

42 

0.0004 

0.0018 

0.0261 

2.5309 

43 

0.0004 

0.0016 

0.0242 

2.2161 

44 

0.0004 

0.0015 

0.0261 

2.4027 

45 

0.0004 

0.0014 

0.0299 

1.9605 

46 

0.0004 

0.0014 

0.0316 

1.9165 

47 

0.0004 

0.0014 

0.0306 

2.2086 

48 

0.0004 

0.0014 

0.0297 

2.3207 

49 

0.0004 

0.0015 

0.0281 

2.1990 

50 

0.0004 

0.0016 

0.0229 

2.8940 

VK  represents  the  diagonal  elements  of  the  ex- 
perimental covariance  matrix  for  the  optimal 
process. 
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Table  IV- 5 


SUBOPTIMAL  VARIANCE  OF  ESTIMATION  ERROR 


K 

SR(1,1) 

SR(2,2I 

VR(1,1I 

VR<2,2> 

1 

0.1667 

2.0000 

0.1308 

1.9041 

2 

0.1018 

1.9833 

0.0899 

1.8888 

3 

0.0777 

1.8799 

0.0725 

1.7924 

4 

0.0682 

1.6698 

0.0615 

1.5970 

5 

0.0653 

1.3827 

0.0627 

1.3288 

6 

0.0653 

1.0666 

0.0578 

1.0320 

T 

0.0661 

0.7655 

0.0603 

0.7471 

8 

0.0662 

0.5098 

0.0610 

0.5028 

9 

0.0650 

0.3136 

0.0603 

0.3134 

10 

0.0622 

0.1780 

0.0598 

0.1806 

11 

0.0579 

0.0954 

0.0553 

0.0990 

12 

0.0525 

0.0538 

0.0532 

0.0549 

13 

0.0464 

0.0406 

0.0490 

0.0393 

14 

0.0399 

0.0443 

0.0406 

0.0407 

15 

0.0336 

0.0558 

0.0351 

0.0506 

16 

0.0277 

0.0687 

0.0294 

0.0627 

17 

0.0225 

0.0791 

0.0241 

0.0729 

18 

0.0180 

0.0852 

0.0187 

0.0793 

19 

0.0143 

0.0864 

0.0149 

0.0813 

20 

0.0114 

0.0P34 

0.0120 

0.0794 

21 

0.0091 

0.0770 

0.0096 

0.0743 

22 

0.0075 

0.0685 

0.0072 

0.0669 

23 

0.0063 

0.0589 

0.0061 

0.0582 

24 

0.0055 

0.0492 

0.0048 

0.0493 

25 

0.0049 

0.0401 

0.0046 

0.0406 

26 

0.0046 

0.0320 

0.0041 

0.0327 

27 

0.0044 

0.0251 

0.0038 

0.0258 

28 

0.0042 

0.0196 

0.0040 

0.0202 

29 

0.0041 

0.0153 

0.0039 

0.0157 

30 

0.0039 

0.0121 

0.0038 

0.0123 

31 

0.0038 

0.0099 

0.0036 

0.0099 

32 

0.0037 

0.0084 

0.0035 

0.0083 

33 

0.0035 

0.0074 

0.0034 

0.0074 

34 

0.0033 

0.0068 

0.0030 

0.0069 

35 

0.0032 

0.0065 

0.0028 

0.0066 

36 

0.0030 

0.0063 

0.0026 

0.0063 

37 

0.0028 

0.0062 

0.0025 

0.0061 

38 

0.0026 

0.0061 

0.0025 

0.0058 

39 

0.0025 

0.0060 

0.0026 

0.0056 

40 

0.0023 

0.0058 

0.0026 

0.0054 

41 

0.0022 

0.0056 

0.0025 

0.0052 

42 

0.0021 

0.0054 

0.0023 

0.0050 

43 

0.0020 

0.0052 

0.0023 

0.0047 

44 

0.0019 

0.0050 

0.0022 

0.0045 

45 

0.0018 

0.0047 

0.0020 

0.0042 

46 

0.0017 

0.0045 

0.0020 

0.0040 

47 

0.0016 

0.0043 

0.0018 

0.0039 

48 

0.0016 

0.0041 

0.0017 

0.0038 

49 

0.0015 

0.0039 

0.0017 

0.0038 

50 

0.0015 

0.0037 

0.0017 

0.0038 

SR  represents  two  of  the  diagonal  elements  of  the  sub- 
opt  iroal  theoretical  covariance  of  estimation  error, 
corresponding  to  the  sub-optimal  estimation  of  x. 


and  x 


2* 


VR  follows  exactly  the  same  pattern  as  SR,  but  it  re- 
fers to  the  experimental  covariance  of  estimation 
error. 
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Table  IV-6 
EXPERIMENTAL  MEANS  OF  ESTIMATION  ERROR 


K 

MEANE1 

MEANE2 

MEANE3 

MEANE4 

MEANE1R 

MEANE2R 

1 

0.0312 

0.0500 

-0.0350 

-0.1543 

0.0312 

0.0500 

2 

0.0263 

0.0436 

-0.0543 

-0.5355 

0.0264 

0.0478 

3 

0.0282 

0.0398 

-0.0893 

-0.7945 

0.0285 

0.0439 

4 

0.0330 

0.0398 

-0.1411 

-0.1811 

-0.8346 

0.0328 

0.0379 

5 

0.0362 

0.0349 

-0.4616 

0.0356 

0.0301 

6 

0.0313 

0.0120 

-0.1482 

-0.1316 

0.0328 

0.0213 

7 

0.0289 

-0.0006 

-0.1229 

0.1844 

0.0317 

0.0120 

8 

0.0041 

-0.0439 

0.1008 

0.1273 

0.0178 

0.0030 

9 

-0.0014 

-0.0429 

0.1368 

0.0809 

0.0142 

-0.0054 

10 

0.0055 

-0.0284 

0.0687 

0.1810 

0.0170 

-0.0124 

11 

0.0150 

-0.0193 

-0.0130 

0.0214 

0.0211 

-0.0176 

12 

0.0031 

-0.0223 

0.0684 

-0.0040 

0.0134 

-0.0213 

13 

-0.0017 

-0.0188 

0.0900 

-0.1130 

0.0093 

-0.0238 

14 

-0.0060 

-0.0134 

0.1014 

-0.2676 

0.0052 

-0.0251 

15 

-0.0075 

-0.0083 

0.0901 

-0.3582 

0.0026 

-0.0254 

16 

-0.0046 

-0.0065 

0.0577 

-0.3345 

0.0031 

-0.0250 

17 

-0.0016 

-0.0069 

0.0383 

0.0395 

0.0042 

-0.0239 

18 

0.0031 

-0.0106 

0.0206 

-0.0433 

0.0073 

-0.0220 

19 

0.0066 

-0.0148 

0.0064 

-0.0677 

0.0101 

-0.0198 

20 

0.0054 

-0.0139 

0.0036 

-0.0764 

0.0085 

-0.0179 

21 

0.0050 

-0.0140 

0.0084 

0.2286 

0.0078 

-0.0164 

22 

0.0024 

-0.0100 

0.0190 

0.1291 

0.0040 

-0.0148 

23 

0.0030 

-0.0111 

0.0238 

0.0811 

0.0053 

-0.0131 

24 

0.0025 

-0.0099 

0.0236 

-0.0649 

0.0046 

-0.0115 

25 

0.0018 

-0.0082 

0.0217 

-0.0134 

0.0033 

-0.0102 

26 

0.0015 

-0.0074 

0.0204 

-0.0375 

0.0030 

-0.0090 

27 

0.0013 

-0.0068 

0.0202 

0.0076 

0.0032 

-0.0080 

28 

0.0011 

-0.0059 

0.0262 

0.1873 

0. 0034 

-0.0069 

29 

0.0006 

-0.0038 

0.0332 

0.1184 

0.0015 

-0.0055 

30 

0.0006 

-0.0025 

0.0299 

-0.2134 

0.0022 

-0.0040 

31 

0.0007 

-0.0016 

0.0225 

-0.1106 

0.0036 

-0.0028 

32 

0.0007 

-0.0008 

0.0181 

-0.0725 

0.0043 

-0.0019 

33 

0.0008 

-0.0001 

0.0093 

-0.2493 

0.0044 

-0.0013 

34 

0.0007 

0.0003 

0.0026 

-0.0454 

0.0033 

-0.0013 

35 

0.0006 

0.0005 

0.0036 

0.0645 

0.0020 

-0.001* 

36 

0.0007 

0.0006 

0.0030 

-0.0675 

0.0020 

-0.0015 

37 

0.0007 

0.0007 

0.0003 

-0.0409 

0.0020 

-0.0016 

38 

0.0007 

0.0006 

-0.0055 

-0.1695 

0.0014 

-0.0019 

39 

0.0007 

0.0002 

-0.0080 

0.0391 

0.0012 

-0.0024 

40 

0.0007 

-0.0001 

-0.0045 

0.0863 

0.0005 

-0.0027 

41 

0.0008 

-0.0002 

0.0006 

0.1190 

0.0015 

-0.0027 

42 

0.0007 

-0.0001 

0.0038 

0.0133 

0.0005 

-0.0023 

43 

0.0008 

0.0002 

0.0048 

0.0356 

0.0015 

-0.0019 

44 

0.0008 

0.0004 

0.0035 

-0.0646 

0.0019 

-0.0014 

45 

0.0008 

0.0005 

0.0012 

-0.0425 

0.0009 

-0.0010 

46 

0.0008 

0.0006 

0.0043 

0.1323 

0.0006 

-0.0007 

47 

0.0007 

0.0009 

0.0081 

0.0145 

-0.0009 

-0.0002 

48 

0.0008 

0.0013 

0.0074 

-0.0337 

-0.0008 

0.0003 

49 

0.0009 

0.0017 

0.0071 

0.0159 

-0.0005 

0.0009 

50 

0.0010 

0.0020 

0.0057 

-0.0529 

0.0001 

0.0015 

MEANE1  through  MEANE4  represent  the  mean  of  estimation 
error  in  the  complete  state  vector  by  means  of  optimal 
filtering. 

MEANE1R  and  MEANE2R  represent  the  sub-optimal  mean  of 
estimation  error  in  states  x-,  and  x9  by  means  of  reduced 
filters.  L  L 
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Fig.  4.3  -  Graphical  plot  of  the  variance  of  estima- 
tion error  in  xi  vs.  time.   Curves  a  and  b  show  the 
experimental  and  theoretical  variance  of  estimation 
error  for  the  optimal  process.   Curves  c  and  d  show 
the  sub-optimal  experimental  and  theoretical  co- 
variance  of  estimation  error. 
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Fig.  4.4  -  Graphical  plot  of  the  variance  of  estima- 
tion error  in  X2  vs.  time.   Curves  a  and  b  show  the 
experimental   and  theoretical  variance  of  estimation 
error  for  the  optimal  process.   Curves  c  and  d  show 
the  sub-optimal  case. 
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Fig.  4.5  -  Graphical  plot  of  the  mean  of  estimation 
error  in  x-^  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  estimation  process 
Curve  b  shows  that  for  the  sub-optimal  case. 
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Fig.  4.6  -  Graphical  plot  of  the  means  of  estimation 
error  in  X2  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  estimation  process. 
Curve  b  shows  that  for  the  sub-optimal  case. 
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V.  IMPLEMENTATION  OF  A  REDUCED  DISCRETE 
KALMAN  FILTER  BY  USING  THE  DOMINANT 
ROOTS  OF  THE  ORIGINAL  PROCESS 

The  general  idea  of  producing  a  reduced  sub-optimal 
filter  (Chapter  III)  or  several  sub-optimal  filters  of  re' 
duced  order  (Chapter  IV),  is  presented  here  using  a  dif- 
ferent approach. 

A.    THE  DOMINANT-ROOTS  APPROACH 

It  is  assumed  as  in  the  previous  chapters,  that  the 
full  system  dynamics  are  represented  by 


x(k+l)  =  $x(k)  +  r  [uM   +   w(k) 


(2.9) 


with  measurement 


z(k)  =  Hx(k)  +  v(k)  , 


(2.3) 


and  with  the  random  processes  involved,  w(k)  and  v(k), 
fully  described  by  Eqs.  (2.11)  through  (2.16). 

It  is  clear  that  the  order  of  the  system  is  n  and 
that  the  goal  of  the  problem  is  to  estimate  only  r  states 
(r<n)  of  the  process;   therefore  it  will  be  necessary  to 
produce  a  reduced  model  of  the  plant  in  order  to  use  its 
dynamics  in  the  implementation  of  the  reduced  filter. 


It  is  common  practice  when  handling  classical  auto- 
matic control  problems  to  assume,  or  to  require,  that  the 
high-order  system  be  dominated  by  one  pair  of  complex  con- 
jugate roots.   This  will  normally  permit  the  use  of  a  se- 
cond-order characteristic  equation,  and  will  simplify 
analysis  and  synthesis  because  the  second-order  system  has 
been  completely  solved  and  all  kind  of  curves,  etc.,  are 
available.   Applying  the  same  idea  to  the  n   -order,  linear 
plant  under  consideration,  but  without  the  constraint  of  a 
pure  second-order  approximation,  the  dominant-roots  repre- 
sentation will  be  determined  mainly  by  the  number  r  of 
states  to  be  estimated  and  by  the  amount  of  performance  de- 
gradation introduced  when  choosing  the  reduced  model. 

It  is  noted  that  the  number  r  of  states  to  be  estimat- 
ed will  specify  the  minimum  number  of  roots  which  can  be 
selected  as  dominant,  but  an  upper  limit  in  the  selection 
of  the  dominant  roots  would  be  determined  by  the  designer's 
judgement,  based  upon  the  amount  of  degradation  introduced 
by  the  reduced  model. 

Due  to  the  qualitative  nature  of  the  root-dominance 
criterion,  it  is  very  difficult  to  state  a  general  proce- 
dure for  the  selection  of  the  reduced-order  model.   There- 
fore a  specific  example  will  be  used  to  illustrate  the 
method . 
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Consider  the  fourth-order  linear  system  with  continir 
ous  dynamics  used  in  Examples  1  and  2  and  described  by 


x(t)  = 


0 

1 

0 

0" 

0 

0 

1 

0 

1 

0 

0 

0 

1 

-400 

-320 

-118 

-19  1 

x(t)  + 


0  l[u(t)+w(t) 

0 
0 


L400] 


with  I.C 


x(0)  =  0  , 


and  measurement 


(5.1) 


(t)  = 


10   0   0 


x 


(t)  +  v(t) 


The  task  of  the  problem  is  to  estimate  only  two  states 
of  the  original  process  (x-,  and  Xo^* 

When  solving  Example  1  (Chapter  3),  the  whole  process 
(Eq.  5.1)  was  discretized  and  then  the  dynamics  of  the  dis- 
crete system  were  partitioned  according  to  the  number  of 
states  to  be  estimated.   Here,  instead  of  discretizing  the 
whole  system,  the  order  of  the  continuous  process  will  be 
lowered  by  proper  choice  of  the  dominant  roots,  and  then 
the  reduced  continuous  system  will  be  discretized. 
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The  system  described  by  Eq .  (5.1)  may  be  rewritten  as 

xx(t)  =  x2(t)  ,  (5.2a) 

x2(t)  =  x3(t)  ,  (5.2b) 

x3(t)  -  x4(t)  ,  (5.2c) 


and 


x4(t)  =  -400x1(t)-320x2(t)-118x3(t)-19x4(t)+400  u(t)+w(t)  . 

(5. 2d) 

For  convenience  a  new  variable  (y(t)  =  x-,(t)  is  defined; 
then  Eqs.  (5.2a,b,c)  become 

(3(t)  =  x2(t)  , 

•  * 

I7(t)  =  X3(t)  , 

£(.t)  =  x4(t)  ,. 

and  by  direct  substitution  of  these  three  equations  into 
Eq .  (5. 2d)  the  latter  becomes 

....        »•»         ■*#         , 

r(t)  +  19  9  (t)  +  118  #(t)  +  320e(t)  +  40(W(t)=40o[u(t)+w(t)j, 

(5.3) 


where 


d49(t) 


dt4 


84 


The  dominant-root  analysis  of  Eq .  (5.3),  which  now  fully 
represents  the  fourth-order  system  under  consideration, 
will  be  made  with  the  help  of  a  pole-zero  representation. 
This  will  require  (as  an  easy  approach)  transformation  of 
Eq .  (5.3)  into  the  s-domain,  accomplished  here  by  use  of 
the  Laplace  transform.   Then,  considering  all  initial  con- 
ditions equal  to  zero,  Eq .  (5.3)  transforms  into 


®(s)s4+19]Ks)s3+118®(s)s2+320e(s)s+400e(s)=400  [u(s)+W(s) 


which  after  proper   factorization   gives 


or 


(ffifO    =  40o[u(s)+W(s)J 

WV J         (s+10)(s+5)(s+2-j2)(s+2+j2)    ' 


®(s) 


400 


U(s)+W(s)         (s+10)(s+5)(s+2-j2)(s+2+j2) 


(5.4) 


Eq.  (5.4),  may  be  interpreted  as  the  transfer  function  of 
the  fourth-order  system  relating  (H)(s)  with  the  forcing 
functions  (inputs). 

A  pole-zero  representation  of  Eq .  (5.4)  (Fig. 5. lb), 
shows  a  pair  of  complex  conjugate  roots  located  at 
s  =  -2±j2,  and  two  real  roots  located  at  s  =  -5  and  s = -10 
respectively.   If  the  dominance  criterion  is  based  upon  the 
size  of  the  real  portion  of  the  roots  under  consideration, 
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a) 


-3S 


b). 


400 


(s+10)(s+5)(sz+4s+8) 


X 


-8     -6 


-4     -2 


X 


-l>®(s) 


Mu 


n 


400 


-2 


c 


uy 


Fig.  5.1  -  a).   Block  diagram  representation  of  the  pro- 
cess described  by  Eq.  (5.4). 

b).  Pole-zero  diagram  for  Eq .  (5.4)  indicating  four  roots 
at  s  =  -5,  s  ■  -10,  s  -  -2±j2  respectively,  and  with  400 
considered  to  be  the  gain  of  the  system. 
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then  it  is  clear  for  this  particular  example,  that  the  less 
dominant  root  will  be  at  s  =  -10  (this  will  give  rise  to  a 
term  of  the  form  £  ),  and  the  most  dominant  roots  will 

be  those  at  s  =  -2+j2,  with  the  root  at  s  =  -5  producing  an 
intermediate  effect. 

As  stated  before,  the  purpose  of  this  problem  is  to 
estimate  only  two  states  of  the  original  system  (x-,  and  X2). 
It  therefore  seems  a  good  approach  to  select  the  roots  at 
s  =  -2ij2  in  order  to  produce  a  reduced  second-order  con- 
tinuous system,  whose  dynamics  will  be  used  in  the  imple- 
mentation of  the  reduced  filter.   It  is  important  to  note 
that  the  gain  of  the  reduced  system  should  be  corrected  in 
order  to  compensate  for  the  disappearance  of  two  of  the 
roots  of  the  original  system.   In  this  case  the  gain  of  the 
reduced  system  was  made  equal  to  the  product  of  the  magni- 
tudes of  the  dominant  roots,  that  is 


1 

gain  reduction  =   -2-j2p-2+j2 

See  Fig.  (5.2) 

Another  important  consideration  should  be  based  upon 
the  fact  that  the  performance  degradation  introduced,  and 
the  number  m  of  inputs  (forcing  functions)  to  the  original 
system,  will  strongly  determine  the  order  of  the  reduced 
model  when  using  dominant  roots. 


). 


b). 


>Or(s) 


2 


-^ 


X 


-2 


Fig.  5.2  -  a).  Block  diagram  representation  of  the 
reduced  model. 

b) .  Pole-zero  diagram  showing  only  the  dominant 
roots.   Gain  of  the  system  has  been  adjusted  to  8. 
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After  the  selection  of  the  pair  of  complex  conjugate 
roots  at  s  =  -2Jlj2,  the  reduced  model  may  be  expressed  as 


0r(s) 


8 


U(s)+W(s) 


(s+2-j2)(s+2+j2)      s2+4s+8 


or  transferring  back  to  the  differential  equation  form, 


Q±t)   +   4^(t)  +  8$ft)  =  8[u(t)+w(t)' 


Recalling  that  u(t)    =   x-,(t),  the  second-order  reduced 
model  becomes 

Xl(t)  =  x2(t) 

x2(t)  =  -8x1(t)  -4x2(t)  +  8[u(t)+w(t)j  , 

or 


xr(t)   = 


0    -1 

-8    -4 


xr(t)   + 


u(t)+w(t)l 


which  represents  the  reduced  model  from  which  the  reduced 
filter  will  be  implemented  (see  Fig.  5.2). 

No  special  derivation  of  the  reduced-filter  equations 
is  necessary  because  the  form  of  the  equations  will  be  ex- 
actly the  same  as  those  for  the  optimal  filtering,  and 
only  the  order  will  change. 
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The  discrete  reduced  model  is  obtained  from  Eq .  (5.5) 


using 


A  T 

$r  -  e  r 


(5.6) 


and 


r 

r 


=  t 


A  T  rT   -AT 
'  r  f     €     r  dTBr  , 
-^  o 


(5.7) 


with  A  (r  x  r)  and  B  (r  x  m)  representing  the  reduced  sys- 
tem and  reduced  distribution  matrices  respectively  for  the 
continuous  case,  and  $  (r  x  r)  and  T  (r  x  m)  being  the 
state  transition  matrix  and  distribution  matrix  for  the 
reduced  discrete  case.   Then  the  discretized  reduced  pro- 
cess can  be  described  as 


xr(k+l)  =  $rxr(k)  +  F  ju(k)+w(k) 
with  measurement 


z(k)  =  Hx(k)  +  v(k)  . 

As  in  the  partitioned  method  presented  in  Chapter  3, 
all  the  measured  outputs  are  considered  to  be  included  in 
the  reduced  model,  and  furthermore  the  H  matrix  is  assumed 
to  be  partitioned  as 


H  = 


Hl  !  H2J 


with  H2  =  0,  and  H..  (p  x  r). 
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Considering  the  problem  of  initialization  for  the  dis- 
crete reduced  covariance  equations,  it  is  solved  exactly  as 
in  Chapter  3,  by  partitioning  the  P(0/-1)  matrix  from  the 
filter,  as 


P(0/-1)  = 


rpn(o/-D  i  p12(o/-d 

p21(o/-i7  :  p227o7-d 


and  choosing 


Pr(0/-1)  -  P-  (0/-1) 


with  P1;L(0/-1)  (r  x  r)  and  P  (0/-1)  (r  x  r) ,  the  latter 
representing  the  initial  value  of  covariance  for  the  re- 
duced filter. 

For  the  reduced  filter  initial  states, 


•N 


£r(o/-i)  =  ^(o/-!) 


where 


x(0/-l)  = 


x^O/.l)' 

x2(0/-l) 
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B.    EXAMPLE  3 

The  same  example  used  in  the  previous  methods  is  em- 
ployed here  to  illustrate  the  method  of  dominant  roots. 

The  same  initialization  values,  same  statistics  for 
the  random  sequences,  and  same  measurements  as  those  used 
in  Example  1,  were  used  in  this  example  for  both  the  full 
and  reduced  systems.   The  corresponding  values  for  the  re- 
duced filter  are  listed  here  for  convenience 


xr(0/-l)  = 


0  I 
0 


Pr(0/-1)  = 


10.5   0 


0    2.0 


and 


h  =  ;i  oj 


v(k)  :  N(0.0,  0.25)  , 
w(k)  :  N(0.0,  0.01)  . 

A  significant  difference  between  the  two  examples  is  the 
determination  of  the  reduced  transition  and  distribution 
matrices  which  in  this  case  were  computed  as  indicated  by 
Eqs.  (5.6)  and  (5.7),  with 


0     1 

i 

1-8    -4 


[o: 


B 
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Due  to  the  similarity  between  the  partitioning  method  and 
the  dominant-roots  approach  for  this  particular  problem, 
the  same  computer  program  was  used  for  the  Monte  Carlo  simula- 
tion and  solution  of  Examples  1  and  3,  with  the  only  change 
introduced  regarding  the  computation  of  the  dynamics  for 
the  reduced  model. 

Table  (V-l)  shows  the  initialization  values  for  the 

reduced  model,  as  well  as  the  values  for  the  $   and  V    ma- 

r      r 

trices.   Tables  (V-2)  through  (V-5)  present  numerical  re- 
sults of  the  simulation,  and  a  graphical  representation  of 
the  simulation  is  shown  in  Figs.  (5.3)  through  (5.6). 

Fig.  (5.7)  shows  a  complete  block  diagram  of  the  re- 
duced filtering  process  for  this  particular  example. 
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Table  V-l 


PHI  MATRIX 

0.9999 

-0.0066 

-0.3644 

-12.2713 

GAMMA  MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 


0.0499 

0.9946 

-0.2981 

-10.1815 


0.0012 

0.0480 

0.8871 

•3.9181 


0.0000 
0.0009 
0.0307 
0.3043 


REDUCED  PHI  MATRIX 


0.9906 
-0.3613 


0.0452 
0.8100 


REDUCED  GAMMA  MATRIX 

0.0093 
0.3613 


Q  MATRIX  (REDUCED) 


0.0000 
0.0000 


0.0000 
0.0000 


VARIANCE  OF  FORCING  NOISE  =0.01 
VARIANCE  OF  MEASUREMENT  NOISE  =0.25 


REDUCED  P(0/-1) 

0.5000     0.0 
0.0        2.0000 

INITIALIZATION  -  X(0/-1) 

0.0 
0.0 
0.0 
0.0 
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Table  V-2 
GAINS    FOR    FULL    AND    REDUCED    FILTERS 


G(  1,1) 


G(2,l) 


G<3, 1) 


G<4,1) 


GR(  1,11 


GR<2,1) 


1 

0.6667 

0.0 

2 

0.4071 

0.2333 

3 

0.3097 

0.4075 

4 

0.2677 

0.5088 

5 

0.2477 

0.5314 

6 

0.2349 

0.4901 

7 

0.2227 

0.4108 

8 

0.2087 

0.3175 

9 

0.1928 

0.2268 

10 

0.1758 

0.1474 

11 

0.  15  8* 

0.0825 

12 

0.1412 

0.0322 

13 

0.1247 

-0.0051 

14 

0.1093 

-0.0315 

15 

0.0950 

-0.0490 

16 

0.0819 

-0.0596 

17 

0.0701 

-0.0647 

18 

0.0596 

-0.0660 

19 

0.0502 

-0.0643 

20 

0.0420 

-0.0606 

21 

0.0348 

-0.0556 

22 

0.0286 

-0.0499 

23 

0.0233 

-0.0438 

24 

0.0188 

-0.0378 

25 

0.0151 

-0.0320 

26 

0.0120 

-0.0266 

27 

0.0095 

-0.0217 

28 

0.0075 

-0.0174 

29 

0.0059 

-0.0136 

30 

0.0047 

-0.0105 

31 

0.0037 

-0.0078 

32 

0.0031 

-0.005  7 

33 

0.0026 

-0.0041 

34 

0.0022 

-0.0028 

35 

0.0020 

-0.0018 

36 

0.0018 

-0.0011 

37 

0.0018 

-0.0006 

38 

0.0017 

-0.0003 

39 

0.0017 

-0.0001 

40 

0.0017 

0.0000 

41 

0.0017 

0.0001 

42 

0.0017 

0.0001 

43 

0.0017 

0.0001 

44 

0.0017 

0.0001 

45 

0.0017 

0.0000 

46 

0.0017 

0.0000 

47 

0.0017 

-0.0000 

48 

0.0017 

-0.0000 

49 

0.0017 

-0.0001 

50 

0.0017 

-0.0001 

0.0 

-0.2095 

•0.6537 

•1.1719 

-1.6080 

•1.8640 

•1.9266 

-1.8378 

-1.6557 

•1.4305 

-1.1970 

•0.9760 

-0.7780 

•0.6072 

•0.4637 

-0.3457 

-0.2505 

-0.1752 

-0.1167 

-0.0723 

-0.0395 

-0.0161 

-0.0001 

0.0100 

0.0157 

0.0181 

0.0181 

0.0166 

0.0141 

0.0111 

0.0080 

0.0050 

0.0023 

-0.0001 

-0.0022 

-0.0038 

-0.0050 

-0.0060 

-0.0066 

-0.0070 

-0.0073 

-0.0074 

-0.0074 

-0.0073 

-0.0072 

-0.0071 

-0.0070 

-0.0068 

-0.0067 

-0.0066 


0.0 

•7.2845 
9.4633 
9.0573 
•7.0288 
4.4015 
1.9563 
•0.0721 
1.1851 
1.9077 
2.2348 
2.2961 
2.1931 
1.9984 
1.7610 
1.5120 
1.2705 
1.0474 
0.8480 
0.6746 
0.5269 
0.4036 
0.3026 
0.2215 
0.1576 
0.1084 
0.0714 
0.0444 
0.0254 
0.0127 
0.0046 
•0.0000 
0.0022 
•0.0028 
•0.0024 
•0.0015 
•0.0004 
0.0007 
0.0017 
0.0024 
0.0030 
0.0033 
0.003* 
0.0033 

coon 

0.0028 
0.0025 
0.0022 
0.0018 
0.0015 


0.6667 
0.4014 
0.2896 
0.2286 
0.1895 
0.1613 
0. 1390 
0.1205 
0. 1045 
0.0903 
0.0778 
0.  0666 
0.0566 
0.0478 
0.0401 
0.0333 
0.0274 
0.0224 
0.0181 
0.0144 
0.0114 
0.0089 
0.0069 
0.0052 
0.0039 
0.0029 
0.0021 
0.0016 
0.0011 
0.0008 
0.0006 
0.0004 
0.0003 
0.0003 
0.0002 
0.0002 
0.0002 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0001 
0.0000 
0.0000 


0.0 

0.0323 
0.0546 
0.0527 
0.0360 
0.0134 
-0.0091 
-0.0286 
-0.0436 
-0.0541 
-0.0605 
-0.0633 
-0.0633 
-0.0612 
-0.0576 
-0.0529 
-0.0477 
-0.0422 
-0.0367 
-0.0314 
-0.0264 
-0.0219 
-0.0179 
-0.0143 
-0.0113 
-0.0088 
-0.0067 
-0.0050 
-0.0037 
-0.0026 
-0.0018 
-0.0013 
-0.0008 
-0.0006 
-0.0004 
-0.0002 
-0.0002 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 
-0.0001 


G  (4x1)  represents  the  matrix  of  gains  for  the  optimal 
filter. 

GR(2xl)  represents  the  matrix  of  gains  for  the  reduced 
filter. 
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Table  V-3 


THEORETICAL  VARIANCE  OF  ESTIMATION  ERROR 


K 

PKC1, 1) 

PK(2,2) 

PK(3,3) 

PK(4,4) 

PRC 1,11 

PR (2, 2) 

1 

0,1667 

2.0000 

2.0000 

4.0000 

0. 1667 

2.0000 

2 

0.1018 

1.9603 

1.7605 

242.6279 

0.1004 

1.3335 

3 

0.07^4 

1.8002 

2.3189 

346.3362 

0.0724 

0.8822 

4 

0,0669 

1.5182 

3.9792 

261.1479 

0.0571 

0.5793 

5 

0.0619 

1.1771 

5.7036 

142.2309 

0.0474 

0.3794 

6 

0.0587 

0.8478 

6.6569 

60.6997 

0.0403 

0.2498 

7 

0.0557 

0.5745 

6.6825 

21.3200 

0.0348 

0.1672 

8 

0.0522 

0.3709 

6.0419 

8.2223 

0.0301 

0.1153 

9 

0.0482 

0.2307 

5.0782 

6.9855 

0.0261 

0.0832 

10 

0.0439 

0.1400 

4.0514 

9.3412 

0.0226 

0.0636 

11 

0.0396 

0.0847 

3.1097 

11.6978 

0.0194 

0.0516 

12 

0.0353 

0.0529 

2.3160 

12.9672 

0.0166 

0.0442 

13 

0.0312 

0.0361 

1.6824 

13.1242 

0.0142 

0.0394 

14 

0.0273 

0.0282 

1.1960 

12.4778 

0.0120 

0.0360 

15 

0.0237 

0.02  55 

0.8337 

11.3650 

0.0100 

0.0332 

16 

0.0205 

0.0254 

0.5710 

10.0516 

0.0083 

0.0308 

17 

0.0175 

0.0264 

0.3855 

8.7171 

0.0069 

0.0285 

18 

0.0149 

0.0275 

0.2581 

7.4696 

0.0056 

0.0261 

19 

0.0126 

0.0283 

0.1736 

6.3652 

0.0045 

0.0237 

20 

0.0105 

0.0286 

0.1199 

5.4260 

0.0036 

0.0212 

21 

0.0087 

0.02  82 

0.0877 

4.6526 

0.0028 

0.0188 

22 

0.0072 

0.0273 

0.0701 

4.0336 

0.0022 

0.0165 

23 

0.0058 

0.0259 

0.0619 

3.5512 

0.0017 

0.0142 

24 

0.0047 

0.0241 

0.0595 

3.1854 

0.0013 

0.0121 

25 

0.0038 

0.0221 

0.0601 

2.9159 

0.0010 

0.0102 

26 

0.0030 

0.0199 

0.0619 

2.7237 

0.0007 

0.0084 

27 

0.0024 

0.0177 

0.0639 

2.5918 

0.0005 

0.0069 

28 

0.0019 

0.0155 

0.0654 

2.5057 

0.0004 

0.0056 

29 

0.0015 

0.0133 

0.0659 

2.4529 

0.0003 

0.0044 

30 

0.0012 

0.0114 

0.0654 

2.4236 

0. 0002 

0.0035 

31 

0.0009 

0.0096 

0.0640 

2.4099 

0.0001 

0.0027 

32 

0.0008 

0.0081 

0.0618 

2.4059 

0.0001 

0.0021 

33 

0.0006 

0.0067 

0.0590 

2.4071 

0.0001 

0.0016 

34 

0.0006 

0.0056 

0.0557 

2.4105 

0.0001 

0.0012 

35 

0.0005 

0.0046 

0.0523 

2.4141 

0.0001 

0.0008 

36 

0.0005 

0.0039 

0.0489 

2.4167 

0.0000 

0.0006 

37 

0.0004 

0.0033 

0.0455 

2.4177 

0.0000 

0.0004 

38 

0.0004 

0.0028 

0.0424 

2.4170 

0.0000 

0.0003 

39 

0.0004 

0.0024 

0.0396 

2.4146 

0.0000 

0.0002 

40 

0.0004 

0.0022 

0.0371 

2.4109 

0.0000 

0.0002 

41 

0.0004 

0.0020 

0.0349 

2.4062 

0.0000 

0.0001 

42 

0.0004 

0.0018 

0.0330 

2.4008 

0.0000 

0.0001 

43 

0.0004 

0.0017 

0.0315 

2.3951 

0.0000 

0.0001 

44 

0.0004 

0.0017 

0.0303 

2.3894 

0.0000 

0.0001 

45 

0.0004 

0.0017 

0.0293 

2.3839 

0.0000 

0.0001 

46 

0.0004 

0.0016 

0.0285 

2.3787 

0.0000 

0.0001 

47 

0.0004 

0.0016 

0.0279 

2.3741 

0.0000 

0.0001 

48 

0.0004 

0.0016 

0.0275 

2.3700 

0.0000 

0.0001 

49 

0.0004 

0.0016 

0.0272 

2.3664 

0.0000 

0.0000 

50 

0.0004 

0.0016 

0.0270 

2.3635 

0.0000 

0.0000 

PK  represents  the  diagonal  elements  of  the  theoretical 

covariance  matrix  for  the  optimal  filter. 

PR  represents  the  diagonal  elements  of  the  theoretical 

covariance  matrix  for  the  reduced  filter. 
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Table  V-4 
EXPERIMENTAL  VARIANCE  OF  ESTIMATION  ERROR 


K 

VK(1,1J 

VK(2t2l 

VK  (  3  , 3  ) 

VK(4,4) 

VR(1, 1) 

VR(2,2) 

1 

0.1510 

1.9066 

1.9733 

4.0880 

0.1486 

1.9041 

2 

0.1022 

1.8632 

1.7310 

225.4341 

0.1008 

1.9624 

3 

0.0807 

1.7518 

2.2373 

326.4941 

0.0801 

2.0096 

4 

0.0667 

1.5112 

3.81.45 

257.3462 

0.0684 

1.9437 

5 

0.0677 

1.1743 

5.6299 

138.0874 

0.0713 

1.7625 

6 

0.0529 

0.8372 

6.0810 

60.2167 

0.0664 

1.5128 

7 

0.0518 

0.5928 

6.5273 

23.2761 

0.0728 

1.2314 

8 

0.0487 

0.3939 

6.020  3 

8.7271 

0.0783 

0.9501 

9 

0.0464 

0.2513 

5.1600 

7.2510 

0.0825 

0.6993 

10 

0.0457 

0.1612 

4.4359 

9.6855 

0.0860 

0.4897 

11 

0.0420 

0.1002 

3.5490 

12.4134 

0.0863 

0.3264 

12 

0.0391 

0.0617 

2.7683 

14.0120 

0.0890 

0.2020 

13 

0.0368 

0.0401 

2.0579 

14.2776 

0.0877 

0.1218 

14 

0.0297 

0.0300 

1.3831 

14.0637 

0.0830 

0.0696 

15 

0.0264 

0.0256 

0.9971 

11.9529 

0.0804 

0.0388 

16 

0.0226 

0.0247 

0.6811 

12.3539 

0.0756 

0.0270 

17 

0.0188 

0.0260 

0.4289 

9.7735 

0.0704 

0.0274 

18 

0.0153 

0.0261 

0.2742 

7.4905 

0.0639 

0.0335 

19 

0.0131 

0.0271 

0.1830 

6.0835 

0.0579 

0.0440 

20 

0.0112 

0.0269 

0.1292 

5.5882 
5.2477 

0.0518 
0.0456 

0.0552 

21 

0.0094 

0.0274 

0.0892 

0.0661 

22 

0.0074 

0.0252 

0.0707 

4.1038 

0.0396 

0.0745 

23 

0.0061 

0.0240 

0.0651 

3.5547 

0.0340 

0.0811 

24 

0.0047 

0.0216 

0.0634 

2.8311 

0.0287 

0.0846 

25 

0.0039 

0.0202 

0.0596 

3.2233 

0.0239 

0.0858 

26 

0.0032 

0.0185 

0.0569 

2.5204 

0.0196 

0.0844 

27 

0.0025 

0.0166 

0.0547 

2.7723 

0.0159 

0.0809 

28 

0.0020 

0.0149 

0.0546 

2.8033 

0.0126 

0.0756 

29 

0.0016 

0.0130 

0.0550 

2.6636 

0.0098 

0.0688 

30 

0.0013 

0.0114 

0.0581 

2.0651 

0.0075 

0.0612 

31 

0.0010 

0.0098 

0.0584 

2.4213 

0.0057 

0.0532 

32 

0.0008 

0.0083 

0.0567 

2. 7746 

0.0042 

0.0453 

33 

0.0007 

0.0071 

0.0539 

2.5252 

0.0031 

0.0379 

34 

0.0006 

0.0060 

0.0490 

2.1855 

0.0022 

0.0313 

35 

0.0005 

0.0050 

0.0435 

2.5705 

0.0016 

0.0254 

36 

0.0005 

0.0042 

0.0420 

2.2536 

0.0011 

0.0204 

37 

0.0004 

0.0035 

0.0437 

2.5342 

0.0008 

0.0162 

38 

0.0004 

0.0029 

0.0446 

2.3070 

0.0006 

0.0127 

39 

0.0004 

0.0025 

0.0421 

2.0512 

0.0005 

0.0098 

40 

0.0004 

0.0022 

0.0366 

2.4148 

0.0005 

0.0075 

41 

0.0004 

0.0019 

0.0300 

2.2901 

0.0004 

0.0056 

42 

0.0004 

0.0018 

0.0256 

2.5146 

0.0004 

0.0043 

43 

0.0004 

0.0017 

0.0241 

2.2052 

0.0004 

0.0033 

44 

0.0004 

0.0016 

0.0261 

2.3978 

0.0004 

0.0025 

45 

0.0004 

0.0015 

0.0300 

1.9640 

0.0005 

0.0020 

46 

0.0004 

0.0014 

0.0317 

1.9171 

0.0005 

0.0017 

47 

0.0004 

0.0014 

0.0307 

2.2013 

0.0005 

0.0015 

48 

0.0004 

0.0014 

0.0299 

2.3214 

0.0005 

0.0015 

49 

0.0004 

0.0015 

0.0281 

2.2047 

0.0005 

0.0015 

50 

0.0004 

0.0016 

0.0229 

2.8952 

0.0005 

0.0016 

VK  represents  the  diagonal  elements  of  the  experimental 
covariance  matrix  for  the  optimal  process. 

VR  represents  the  diagonal  elements  of  the  experimental 
covariance  matrix  for  the  reduced  process. 
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Table  V-5 


EXPERIMENTAL  MEANS  OF  ESTIMATION  ERROR 


K 

MEANE1 

MEANE2 

MEANE3 

MEANE4 

MEANE1R 

MEANE2R 

1 

0,0496 

0.0500 

-0.0350 

-0.1543 

0.0496 

0.0500 

2 

0.0372 

0.0392 

-0.0572 

-0.6276 

0.0349 

-0.1157 

3 

0,0356 

0.0307 

-0.0903 

-0.7992 

0.0266 

-0.2315 

4 

0.0380 

0.0271 

-0.1340 

-0.7676 

0.0193 

-0.3045 

5 

0.0395 

0.0203 

-0.1637 

-0.3713 

0.0100 

-0.3441 

6 

0.0333 

-0.0030 

-0.1219 

-0.0534 

-0.0047 

-0.3578 

7 

0.02*59 

-0.0146 

-0.0909 

0.2345 

-0.0167 

-0.3498 

8 

0.0044 

-0.0565 

0.1352 

0.1478 

-0.0403 

-0.3249 

9 

-0.0017 

-0.0536 

0.1711 

0.0767 

-0.0515 

-0.2894 

10 

0.0049 

-0.0373 

0.1013 

0.1589 

-0.0548 

-0.2527 

11 

0.0142 

-0.0265 

0.0168 

-0.0124 

-0.0556 

-0.2173 

12 

0.0021 

-0.0281 

0.0951 

-0.0444 

-0.0653 

-0.1749 

13 

-0.0028 

-0.0232 

0.1136 

-0.1564 

-0.0704 

-0.1356 

14 

-0.0072 

-0.0168 

0.1219 

-0.3115 

-0.0742 

-0.0982 

15 

-0.0086 

-0.0108 

0.1078 

-0.4010 

-0.0757 

-0.0651 

16 

-0.0058 

-0.0083 

0.0729 

-0.3752 

-0.0745 

-0.0382 

17 

-0.0028 

-0.0080 

0.0511 

0.0015 

-0.0724 

-0.0154 

18 

0.0020 

-0.0112 

0.0313 

-0.0783 

-0.0691 

0.0028 

19 

0.0055 

-0.0150 

0.0153 

-0.0996 

-0.0657 

0.0183 

20 

0.0044 

-0.0137 

0.010Q 

-0.1052 

-0.0635 

0.0337 

21 

0.0040 

-0.0136 

0.0142 

0.2029 

-0.0606 

0.0456 

22 

0.0014 

-0.0094 

0.0236 

0.1062 

-0.0581 

0.0570 

23 

0.0022 

-0.0103 

0.0273 

0.0610 

-0.0543 

0.0640 

24 

0.0017 

-0.0089 

0.0262 

-0.0824 

-0.0506 

0.0699 

25 

0.0011 

-0.0072 

0.0235 

-0.0285 

-0.0469 

0.0738 

26 

0.0008 

-0.0063 

0.0215 

-0.0504 

-0.0430 

0.0758 

27 

0.0007 

-0.0056 

0.0207 

-0.0033 

-0.0390 

0.0762 

28 

0.0006 

-0.0048 

0.0263 

0.1783 

-0.0351 

0.0755 

29 

0.0001 

-0.0027 

0.0329 

0.1111 

-0.0314 

0.0744 

30 

0.0002 

-0.0014 

0.0292 

-0.2192 

-0.027^ 

0.0723 

31 

0.0003 

-0.0006 

0.0215 

-0.1151 

-0.0241 

0.0693 

32 

0.0004 

0.0002 

0.0170 

-0.0758 

-0.020^ 

0.0656 

33 

0.0005 

0.0008 

0.0080 

-0.2516 

-0.0175 

0.0613 

34 

0.0004 

0.0011 

0.0012 

-0.0467 

-0.0146 

0.0564 

35 

0.0004 

0.0012 

0.0022 

0.0639 

-0.0119 

0.0513 

36 

0.0005 

0.0014 

0.0016 

-0.0675 

-0.0094 

0.0463 

37 

0.0006 

0.0014 

-0.0011 

-0.0404 

-0.0072 

0.0413 

38 

0.0006 

0.0012 

-0.0069 

-0.1685 

-0.0053 

0.0363 

39 

0.0006 

0.0007 

-0.0093 

0.0404 

-0.0036 

0.0313 

40 

0.0006 

0.0003 

-0.0057 

0.0877 

-0.002? 

0.0266 

41 

0.0007 

0.0002 

-0.0006 

0.1207 

-0.0009 

0.0225 

42 

0.0007 

0.0003 

0.0027 

0.0151 

0.0001 

0.0139 

43 

0.0008 

0.0004 

0.0038 

0.0374 

0.0010 

0.0157 

44 

0.0008 

0.0006 

0.0026 

-0.0628 

0.0017 

0.0129 

45 

0.0008 

0.0007 

0.0004 

-0.0407 

0.0022 

0.0103 

46 

0.0008 

0.0007 

0.0035 

0.1341 

0.0027 

0.0080 

47 

0.0007 

0.0010 

0.0075 

0.0162 

0.0031 

0.0062 

48 

0.0008 

0.0014 

0.006*5 

-0.0321 

0.0033 

0.0048 

49 

0.0009 

0.0017 

0.0066 

0.0174 

0.0035 

0.0036 

50 

0.0010 

0.0020 

0.0053 

-0.0515 

0.0037 

0.0026 

MEANE1  through  MEANE4  represent  the  mean  of  estimation 
error  in  the  complete  state  vector  by  means  of  optimal 
filtering. 

MEANE1R  and  MEANE2R  represent  the  mean  of  estimation 
error  in  states  x,  and  Xj   by  means  of  a  reduced  filter 
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0.35f 


Time  in  seconds 

Fig.  5.3  -  Graphical  plot  of  the  variance  of  estima- 
tion error  in  xi  vs.  time.   Curves  a  and  b  show  the 
theoretical  variance  of  estimation  error  for  the 
reduced  and  full  filters  respectively.   Curves  c  and 
d  show  the  experimental  variance  for  the  full  and 
reduced  processes  respectively. 
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Time  in  seconds 

Fig.  5.4  -  Graphical  plot  of  the  variance  of  estima- 
tion error  in  xo  vs.  time.   Curves  a  and  b  show  the 
theoretical  variance  of  estimation  error  for  the 
reduced  and  full  filters  respectively.   Curves  c  and 
d  show  the  experimental  variance  for  the  full  and 
reduced  processes  respectively. 
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Fig.  5.5  -  Graphical  plot  of  the  means  of  estimation 
error  in  x;l  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  Kalman  filter,  and 
curve  b  shows  that  for  the  reduced  filter. 
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Fig.  5.6  -  Graphical  plot  of  the  means  of  estimation 
error  in  X2  vs.  time.   Curve  a  shows  the  mean  of 
estimation  error  for  the  optimal  Kalman  filter,  and 
curve  b  shows  that  for  the  reduced  filter. 
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VI.  DISCUSSION  OF  RESULTS 

The  three  different  approaches  presented  in  this 
thesis  have  as  a  common  objective  a  decrease  in  the  compu- 
tational requirements  involved  in  handling  systems  of  high 
dimension,  based  upon  the  fact  that  only  a  preselected  re- 
duced number  (r)  of  states  is  of  significant  importance  in 
the  estimation  process. 

For  the  particular  example  chosen  as  an  illustration 
for  the  three  methods,  the  partitioning  idea  discussed  in 
Chapter  III,  although  it  is  the  easiest  to  apply,  appears 
to  give  the  least  attractive  results  in  the  reduced  esti- 
mation of  states  X-,  and  x^   where  a  small  bias  is  produced. 
This  seems  to  be  a  direct  consequence  of  the  degradation 
introduced  when  $  -.  (which  in  this  example  has  relative 
large  elements  as  compared  with  $•,  )  was  not  considered  to 
be  a  part  of  the  dynamics  of  the  primary  system. 

In  the  pseudo-inversion  approach,  implementation  of 
four  scalar  filters,  plus  the  selection  of  a  proper  set  of 
A.  matrices  and  its  pseudo-inverses  A.  was  necessary,  but 
since  it  takes  less  computing  capacity  to  manipulate  four 
scalar  filters  than  one  (4  x  4)  filter,  the  estimation 
process  was  clearly  simplified. 
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The  selection  of  the  partitioning  matrices  A^  is  arbi< 
trary  for  any  particular  problem,  therefore  for  this  ex- 
ample another  selection  was  made  as 


Ax  =  J  1    1   -1   -1  , 


and 


with 


A2  =  ,1 


A3  =  ,1   -1   -1 


I  . 


A,  =  |  1 


4  = 


r 


0.25 


0.25' 
:-0.25: 


L 


-0.25 


A+  - 


0.25  | 
0.25  | 
0.25   ' 
0.25j 


a;  = 


ro. 

25] 

|-0. 

25| 

!-0, 

25 

L°- 

25 

A 


r  0.25] 
+  -0. 25  J 
4  ,  0. 25 | 

1-0.25' 

L 


J 


This  particular  choice,  however,  leads  to  four  fil- 
ters for  which  the  sub-optimal  variance  of  estimation  er- 
ror, both  theoretical  and  experimental,  increases  linearly 
with  time,  which  is  undesireable .   Therefore  this  parti- 
cular partitioning  was  not  used. 

The  plant  for  all  three  examples  was  selected  with  a 
pair  of  clearly  dominant  complex  roots  in  order  to  obtain 
a  good  second-order  approximation  of  the  complete  process 
from  which  the  reduced  filter  could  be  implemented  using 
the  method  presented  in  Chapter  V. 
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APPENDIX  A 

MATRIX  PSEUDO  INVERSION 

A  review  of  matrix  pseudo-inversion  is  presented  in 
this  appendix  as  an  extension  to  the  idea  of  pseudo-inver- 
sion introduced  in  Chapter  IV. 

A  matrix  A   is  said  to  be  the  pseudo-inverse  of  a  rec- 
tangular matrix  A  if 


A  A+A  =  A 


(A-l) 


where  if  A  is  a  square  non-singular  matrix,  i.e.,  A~   ex- 
ists, then 


A  A_1A  =  A 


is  a  trivial  solution  since  A  A    =  I. 

Another  important  definition  may  be  introduced  using 
a  set  of  four  matrix  equations 


A  A+  A  =  A 


A+A    A+        A+ 

AAA   =  A 


A  A 


+ 


T  -  AA+, 


(A-2a) 
(A-2b) 
(A-2c) 


and 


fA+A 


T      + 
=  A  A 


(A-2d) 
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This  definition  implies  that  the  inverse  of  a  matrix  A  is 
established  as  the  unique  pseudo-inverse  if  and  only  if  it 
satisfies  the  set  of  Eqs.  (A-2).   Furthermore,  as  a  conse- 
quence of  the  latter  definition,  the  pseudo-inverse  among 
other  properties  has  the  following: 


and 


A++  =  A 

AT+  =  A+T 

A   =  A   ,  when  A  is  a  non-singular  matrix, 


0+  =  0T 


The  pseudo-inverse  A   of  a  matrix  A  may  be  shown  to 


be  unique  as  follows. 

Let  A.  be  a  set  of  pseudo-inverse  of  A.  where 
i  =  1,2  . . .  M,  and  assume  that  the  A.  are  chosen  so  that 


M 

^A+A.  =  I  .  (a-3) 

i  =  l 

Further  assume  that  another  set  of  pseudo-inverse  Y.  4   A. 

r  1  '       i 

exists  for  each  A^,  and  that  this  new  set  of  Y.  also  sat- 
isfies 


M 

ZYtAi  = x  •  <a-4> 

i=l 
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Then  it  follows  by  proper  use  of  Eqs.  (A-3)  and  (A-4)  that 


M 
i=l 


Y+)A. 
i7  l 


=  0  , 


and  since  the  choice  of  A.  is  arbitrary,  it  follows  that 


Y+  =  A+ 

1         1 


with   i  =  1,2  . . . ,  M 


and  the  uniqueness  property  is  established. 

For  more  complete  treatment  of  this  subject  see 
Penrose  (Ref.8),  Werther  (Ref.9)  and  Deutsch  (Ref.10). 
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APPENDIX  B 

C      IMPLEMENTATION  OF  A  REDUCED  SUBCPTIMAL  LINEAR  FILTER 

C      FCR  DISCRETE  DYNAMICS  PROCESS,  BY  PARTITIONING  THE 

C      DYNAMICS  CF  THE  ORIGINAL  PROCESS.  THE  DEGRADATION  ERR 

C      OP  OF  THE  FILTER  IS  COMPUTED  AND  USED  AS  A 

C      MEANS  OF  CORRECTING  THE  FILTER  VARIANCES.  A  MONTE  CAR 

C      LC  SIMULATION  OF  THE  PRCCESS  IS  PERFORMED  IN  CRDER  TO 

C      OBTAIN  APROXIMATED  VALUES  FOR  THE  MEANS  AND  COVARIANCE 

C      CF  ESTIMATION  ERROR. 

C 

C 

CI  PENSION  GAM(4,1),GAMR(2,1),PHI(4,4).PHIR(2,2),A(4,4) 
l,H(l,4).B(4,l),HT(4,l),AT(4f4),CR(2.2),VW(l,I),Q(4,4) 
2,PK(4,4),HR(l,2),HTR(2,l),G(4,I),U(4,4l ,PK1R(2.2) 
3,PKl(4,4),PKR(2.2),VWR(l,l),GR(2,l),PKOR(2,2),COR(2,2) 
4,E(2,2) .AX  1(4, 4) ,AX2(4,4I ,AX3(4,4) ,BX1 (4, 1 ) , BX2 (4, 1 ) 
5,BX3(4, 1),GAMU(4,1» ,DX2( 1,41, DEG( 2, 21,0X3(2,2) ,UR( 2, 2) 
6,CX4(2,2),A1(2,2I ,A2( 2,2  I ,0X1 (2,1), G AMUR (2,1 »,CX 3(2,1) 
7,CX1(1,1),CX2( 1,1),R(1,1) ,Y1(50,4) .  Y2  (  50  ,4  )  ,  Y3  (50  ,2  ) 
8  ,  Y4  (  50  ,  2  ) ,  Y 5  (  50  ,  2  ) ,  Y6  (  50 ,  4  )  ,  Y7 (  50 ,  2  ) ,  Y8  <  50  , 4  )  ,  Y9 (  50  ,  2 ) 
?,Y10(50,2) ,Y11(50,1),Y12(50,1),XE(4),XER(2),SE1(50>, 
1SE2(5C),SE3(50)1SE4(50),SE1R(50)ISE2R(50)  .SV1K50). 
2SV22(50),SV33(50). SV44( 50 ) , SVR1 1 1  50 ) ,SVR22< 51 ) , X( 5* ) , 
3LL(1),MM(1),LLL(1) ,MMM( 1 ) ,PD1 ( 5C ) ,PD2( 50 ) , VD1 (50  ) 
4 , V  D2 ( 50 ) , PR  1 ( 50 ) , P  R  2 ( 50 ) , VR 1 ( 50 ) , VR  2 ( 50 ) , E  D 1 ( 50 ) , 
5EC2(50) ,ER1(50),ER2(50) 

REAL*8  IT( 12)/»VARIANCES  OF  ESTIMATION  ERROR  FOR  X(l) 
1  LARA  BCX11  PART.  CORR.  E=P  •/ 

REAL*8  IP( 12)/«VARIANCES  OF  ESTIMATION  ERROR  FOR  X(2) 
1  LARA  BCX11  PART.  CORR.  E=P  •/ 

REAL*8  IE(12)/»MEAN  OF  ESTIMATION  ERROR  FOR  X(l) 
1  LARA  B0X11  PART.  CORR.  E=P  •/ 

REAL*8  IF(12)/«MEAN  OF  ESTIMATION  ERROR  FOR  X(2) 
1  LARA  B0X11  PART.  CORR.  E=P  •/ 

REAL  LAB/4H     / 
CAT A  NT,NS,NE,M,MR,N,NR,NZ, I  X, T/40, 50 ,2C0 , 1, 1,4, 2, 1 , 3, 

DATA  A(1!1),A(1.2),A(1,3),A(1,4),A(2,1),A(2,2),A(2,3), 
2M2,4),A(3,1),A(3,2),A(3,3),A(3,4),A(4,1),A(4,2).A(4,3 
3) ,A(4.4)/0.,1. ,0.,0.,0.,0.,I.,0.,0.,0.,0.,l.,-400.,-32 
40. ,-118. ,-19./ 

DATA  B(l,l ) ,8(2,1 ),B( 3,1), B( 4,1 )/0.,0.,0., 400./ 

DATA  H( 1,1)  ,H(1,2) ,H(1,3)  ,H( 1 ,4 )/l . ,0. ,0. ,0./ 

DATA  HT(1,1),HT(2, 1 ) ,HT( 3, 1 ) ,HT( 4, 1 ) /l. ,0. ,0. ,0. / 

DATA  VW(1,1),RR,VWR( 1,1), UK/0. 01, 0.25. 0.01,0. 5/ 

DATA  SVl,Shl,AWl,KADA,LC0R/0.5,0.1,3.O,l,5/ 

CAT A  SI, S2, S3, S4/0. 707, 1.41, 1.41,2.0/ 

DATA  AM1,AM2,AM3,AM4/0.0,0.0,0.0,0.0/ 

WRITE(6,222) 

CC  31  K=1,NS 

SV33(K)=0.C 

SV44(K)=C.C 

SVR22(K)=C.Q 

SE3(K)=0.0 

SE4(K)=C0 

SE2R(K)=0.C 

SV11(K)=0.C 

SV22(K)=0.C 

SVR11(K)=C0 

SE1(K)=C.C 

SE2(K)=C.C 

SE1R(K)=0.C 
31    CONTINUE 
C 
C      MONTE  CARLC  SIMLLATICN 


C 


CC  1  L=1,NE 

CALL  GAUSS(IX,S1,AM1,X1) 

CALL  GAUSS(IX,S2,AM2,X2) 

CALL  GAUSS(IX,S3,AM3,X3) 

CALL  GAUSS(IX,S4,AM4,X4) 

X1K=X1 

X2K=X2 
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49 


X3K=X3 
X4K=X4 

GENERATING  PHI  ANO  GAMMA 


66 
4T 


81 
6 

123 


124 


IFU.GT 

CALL  UN 
CC  49  I 
CC  49  J 
LP( I ,J> 
CALL  UN 
CALL  UN 
CALL  SM 
CALL  SM 
CALL  SM 
CC  66  I 
XL=I 
RP=1.0/ 
PP1=1.( 
CALL  GM 
CALL  SM 
CALL 
CALL 
CALL 
CALL 
CCNTINU 
WRITE(6 
CC  AC  I 
fcPITE<6 
WPITE(6 
WPITE<6 


GM 
GM 

SM 

GM 


•  l)GO  TC  6 

IT(L,N) 

=  1,NR 

=  1,NR 

=  U( I, J) 

ITf PHItNl 

IT( AX1,N) 

PY<  B,T,GAM,N,M,S) 

PY(B,T,BX2,N,M,0) 

PY( A,T,AT,N,N,0) 

=  1,NT 

XL 
/IXL+l.O) 

PPD( AX1,AT,A,N,N,N 

PY( A,RP, AX1 ,N,N,G) 

ADD(PHI,AX1,PHI,N,N) 

PRD(AT,BX2,B,N,N,M) 

PY(8,RP1,BX2,N,M,0) 

ADD(GAM,EX2,GAM,N,M) 

E 

,13C) 

=  1.N 

,125)  (PHIU  ,J)  ,J=1,N) 

,131) 

,125)GAM 


) 


GENERATING  Q 


H2 


CALL    GM 

CALL    GM 

CALL    GM 

V*PITE(6 

CC    81    I 

KPITE(6 

WRITEC6 

bRITE<6 

CCNTINU 

CC    123 

CC    123 

PK1(I,J 

CCNTINU 

PK1(3,3 

PK1(4,4 

PKU  1,1 

PKK2.2 

XE(1)=C 

XE(2)=0 

XE(3)=0 

XE(4)=C 

XER(  1)  = 

XEP(2)= 

CC    124 

00    124 

PK1RU, 

E(I,J)= 

CCNTINU 

PK1R(1, 

PK1R<2, 

E(l,l)= 

E<2,2)= 

IF(L.GT 

*RITE(6 

CC    142 

V«RITE(6 

WRITE(6 

V<PITE(6 


PRO 
TRA 
PRD 
,13 

:h 

,13 

i=i 
j=i 

)=C 

E 

)=2 

)=4 

)=0 

)=2 

• 

.0 

.0 

.0 

o.c 
o.c 
1=1 
J=l 
J)= 

CO 
E 

1)  = 

2)  = 
0.5 
2.0 
•  1) 

\H 

,12 
,13 
,12 


(GAM,VW,BX1,N,M,M) 

(GAM,DX2,N,M) 

(BX1,CX2,Q,N,M,N) 

2) 

N 

5)(Q( I,J),J=1,N) 

5) 

6) 

iN 
,N 
.0 


.5 


,NR 

,NR 

e.o 


0.5 

2.0 


GO    TO    51 

4) 

,N 

5MPKHI  ,J)  ,J  =  1,N) 

5)XE 


110 


143 
51 
C 


12 
13 

14 


WPITEC6.138) 

CC    143    1=1, NR 

WPITE(6,12  5)(E(I,J),J=1,NR) 

CCMINUE 

CC    3    K=1,NS 

CALL    GAUSS(IX,SV1,0.G,V1) 

CALL    GAUSS(IX,SH1,AW1,W1) 

IF(K.EQ.l)GO    TO    4 

X1K=PHI(1, 1)*X1*PHI(1,2)*X2*PHI (1,3>*X3+PHI(1,4)*X4 
1+GAM<1,1)*<UK+W1) 

X2K=PHI<2, 1>*X1+PHI(2,2)*X2+PHI (2,3I*X3*PHI(2,4)*X4 
2-i-GAM<2,l>*<UK*Wl) 

X3K=PHI  <3,1)*X1+PHI  <3,2)*X2-»-PHI  (3  ,  3  )*X3+PH  I  (  3  ,4  )  *X4 
3+GAM(3.1>*(UK+Wl) 

X4K=PHI<4,  1)*X1+PHI<4,  2)*X2*PHI  (  4,  3  )  *X3*PH  I  (  4  ,  4  )  *X4 
4+GAM<4,l>*(UK+*l) 

MEASUREMENT 

Z=Vi  1,1)*X1K+V1 

GAIN  AND  CCVARIANCES  CF  EST,  ERROR  FOR  COMPLETE  SYSTEM 

C/SLL  GMPRD(PK1,FT,BX1,N,N,NZ) 

CALL  GMPRD(H,BX1,CX2,NZ,N,NZ) 
SC=1.C/(CX2(1,1)+RR) 

CALL  SMPY(EX1,SC,G,N,NZ,0 ) 

CALL  GMPRD(G,H,AX1,N,NZ,N) 

CALL  GMSUB(U,AX1,AX1,N,N) 

CALL  GMPPD(AX1,PK1,PK,N,N,N> 

CALL  GMTRACPHI  ,AX2,N,N) 

CALL  GMPRD(PK,AX2,AX1,N,N,N) 

CALL  GMPRDtPHI ,  AX1 f AX2 1 Nt N»N1 

CALL  GMADD<AX2,C,PK1,N,M 

LSING  KALMAN  FILTER  ECUATION  WITH  FULL  PLANT 


IF(K.GT 
CALL  GM 
GO  TO  8 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

E1=X1K- 

E2=X2K- 

E3=X3K- 

E^=X4K- 


GM 
SM 
GM 
GM 
GM 
GM 
SM 
GM 


.l)GO  TO  7 

PRD(U,XE,BX1  ,N,N,M) 

PRO(PHI,XE,BXl,N,N,M) 

PY(GAM,UK,GAMU,N,M,0) 

A0D(BX1,GAMU,BX1,N,M) 

PRD(G,H,AX1,N,M,N) 

SUB(U,AX1,AX2,N,M 

PRD(AX2,eXl,BX2,N,N,M) 

PY(G,Z.BX1,N,M,5) 

ADD(BX1,BX2,XE,N,1) 

XE(  1) 

XE(  3) 
XE<4> 


GAIN  AND  CCVAR.  OF  ESTIMATION  ERROR  FOR  REDUCED  SYSTEM 

IFU.GT.DGC    TC    11 
CC    12    1=1, NR 
CC    12    J=1,NR 
PHIRd  .  J)  =  PHI(  I, J) 
CO    13    1=1, NR 
CC    13    J=1,MR 
GAMR(I,J)=GAM( I, J) 
CC    14    1=1, NR 
CC    14    J=1,NZ 
HR(J,I )=H< J, I) 
FTR(  I, J)  =  HT( I, J) 

C    RECUCED 

CALL    GMPRD(GAMR,VWR,DX1,NR,MR,MR) 


111 


44 
11 


1  1  )r' 
41 


1J 


21 
2000 


37 

17 


CALL  GMTRA 
CALL  GMPPD 
CC  44  1=1, 
CC  44  J=l, 
PKCRd  ,  J)  = 
CALL  GMPRO 
CALL  GMPRO 
CX2< 1,11=1 
CALL  GMPRO 
CALL  GMPRO 
CALL  GMSUB 
IF(LCOR.EC 
IF<K.GT.1> 
CALL  GMPRO 
CONTINUE 
CALL  GMPPD 
CALL  GMSUB 
CALL  GMTRA 
IF(LCOR.EQ 
IF(K.EQ.l) 
CALL  GMPPD 
CALL  GMPRO 
CALL  GMPRO 
DC  20  1=1, 
CC  2C  J=l, 
PKR(  I, J  >  =  P 
E(I, J)=CCR 
CCMINUE 
CALL  GMPRO 
CALL  GMPRO 
CALL  GMADD 


(GAMR,CX2,NR,MR) 

(DX1,CX2,CR,NR,MR,NR) 

NR 


GO  TO  41 
<DX3,E,CCR,NR,NR,NR) 


(DX3,PK1R,PKP,NR,NR,NR) 

<FKR,PKaR,DEG,NR,NR) 

(PHIR,DX4,NR,NR) 

.DGC  TO  20^0 

GO  TO  10 

(PHIR,E,A1,NR,NR,NR) 

( A1,DX4,A2,NR,NR,NR) 

(DX3,A2,CCR,NR,NR,NR) 

NR 

NR 

KR( I, JI+CCRC  I, J) 

(I, J) 

<PKR,CX4,CX3,NR,NR,NR) 

(PHIR,DX3,DX4,NR,NR,NR) 

(DX4,CR,PK1R,NR,NR) 


LSING  KALMAN  FILTER  EQUATION  WITH  RECUCED  PLANT 


IF(K 
CALL 
GO  T 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
E1R= 
E2R  = 
X1  =  X 
X2  =  X 
X3  =  X 
X4=X 


.GT 
GM 

0  1 
GM 
SM 
GM 
GM 
GM 
GM 
SM 
GM 

X1K 

X2K 

IK 

2K 

3K 

4K 


•  DGO    TO    37 

PRD(UR,XER,CX3,NR,NR,1) 

7 

PRD(PHIR,XER,CX1,NR,NR, 1) 

PY(GAMR,UK,GAMUR,NR,M,0  ) 

ADD ( DX1,G AMUR, CX3,NR,1) 

PRD<GP,HP,DX3,NR,1,NR> 

SUB(UR,DX3,CX4,NR,NR) 

PRD(DX4,CX3,CX1,NR,NR,  1) 

PY(GR,Z,CX3,NR,NZ,0) 

ADD(CX3,DX1  ,XER,NR,1) 

-XER(l) 

-XER(2) 


ENSEMBLE  MEAN  OF  PK/K  AND  PK/K  REDUCED,  ANC  GAINS. 


IF(L.GT 
Y1(K,1) 
Y2(K,1) 
Y3<K,1) 
Y1(K,3) 
Y8(K,1) 
Y8(K,2) 
Y8(K,3) 
Y8(K,4) 
Y2(K,3) 
Y«5(K,1) 
Y1C(K,1 
Y1C(K,2 


.DGO  TO 
=  PKU,1) 
=PK(2,2) 
=PK<3,3) 
=PKR(1,1 
=G( 1,1) 
=  G(2, 
=G(3, 
=  G(4, 

=PKP(2,2 
=PK(4,4) 
)=GR(1, 1) 
)=GP(2,1) 


26 


1) 
1) 
1) 


26 


RE*L  VARIANCE  OF  EST.  ERROR  FOR  FULL  AND  REDUCED  SYST. 

SV11(K)=SV11(K)+E1**2 
SV22(K)=SV22(Kl4E2**2 


112 


42 


27 


88 
94 

7 
1 


1111 


77 


SV33(K)=SV33(K)+E3**2 
SV44(K)=SV44(K)4E4**2 
SVR11(K)=SVR11 (K)*E1R**2 
SVR2  2(K)=SVR22(KH-E2R**2 

PEAN  OF  ESTIM.  ERRCR  FCR  FULL  AND  REDUCED  PLANT 

SE1IK)-SE1(K)+E1 

SE2(K)=SE2(KUE2 
SE3(K)=SE3(K)+E3 
SE4(K)  =  SE4(KH-E4 
SE1R<K)=SE1R(K)+E1R 
SE2R(K)=SE2R(K)4E2R 


SYSTEM    DEGRADATION 


3000 


2) 
1) 
2) 

TC    94 


IFUCOR.EQ.l  )GC    TO 

IF(L.GT.1)G0  TO  27 

Y6(K,1J=DEG(1, 1) 

Y6(K,2)=CEG(2 

Y6(K,3)=C0R( 1 

Y6(K,4)=C0R(2t 

CONTINUE 

IF1L.LT.NE  )G0 

DC  88  K1=1»NS 

Y1(K1,2)=SV1UK1)/NE 

Y2(K1,2)=SV22(K1»/NE 

Y3(K1,2)=SV33(K1 )/NE 

Y9(K1,2>=SV44<K1)/NE 

Y4(K1,1)=SEKK1)/NE 

Y5(K1,1)=SE2(K1 )/NE 

YMK1,2)=SE1R(K1>/NE 

Y5(K1,2)=SE2R(K1)/NE 

Y1(K1,A)=SVR11(K1)/NE-Y4(K1,2)**2 

Y2(K1,4)=SVR22(K1)/NE-Y5(K1,2)**2 

Y1HK1,1)  =  SE3(K1)/NE 

Y12(K1,1)  =  SEM  K1J/NE 

X(K1)=K1 

CONTINUE 

CONTINUE 

CONTINUE 

DC    1111    K=l 

PD1(K)=YHK 

PP1(K)=Y1(K 

VDltK)-Yl<  K 

VR1(K)=Y1(K,4) 

PC2(K)=Y2(K,1) 

PR2(K)=Y2(K,3) 

VC2(K)=Y2(K,2) 

VR2(K)=Y2(K,4) 

ED1(K)=YA(K,D 

EC2(K)=Y5(K,1) 

ER1(K)=Y4(K,2) 

EP2(K)=Y5(K,2) 

CONTINUE 

V«RITE<6,222) 

WPITE(6,  16*"  ) 

V«PITE<6,555) 

CO    3CC    K  =  l  ,NS,KADA 

V*PITE(6,lll)K,Y8(K,l),Y8(K,2),Y8(K,3)tY8(K,4),Ylj(K,l) 
ltYlC (K,2) 

CONTINUE 

U RITE (6, 22 2) 

WRITE(6,161) 

WRITE(6,666) 

CO    77    K=1,NS,KACA 

URITE(6,lll)K,YHK,n  ,Y2(K,1)  ,Y3(K,1)  tY9(K,  1  )  fYl(K,3» 
1  •  Y2(  Kf  3  ) 

CCNTINUE 

WRITE(6»222) 

URITE<6, 162) 

WFITE(6,777) 


,  NS 
,1) 
,3) 
♦  2) 
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CC    888    K=1,NS,KA0A 

WRITE<6,111)K,Y1(K,2),Y2<K,2),Y3<K,2>  ,Y9(K,2)  ,Y1(K,A) 
1,Y2<K,A) 

688         CONTINUE 

WRITE(6,222) 
fcRITE<6,163) 
WRITE<6,999) 
CC    757    K  =  l  ,NS,KADA 

V«RITE<6,111)K,YA(K,1)  ,  Y5  (  K  ,  1 )  ,  Y 11  (  K  ♦  1 )  ,  Yl  2  <K  ,  1 )  ,  YA<  K,  2 
1),Y5(K,2> 

757    CONTINUE 

IF(LCOR.EQ.l)GC  TO  4000 

WPITE<6,222) 

fePITE(6,16A) 

WRlTE(6,99e) 

CC    778    K=1,NS,KADA 

bPITE(6,lll)K,Y6(K,l),Y6<K,2),Y6<i',3),Y6<K,A) 

778         CONTINUE 

4000      CCMINUE 

IF<LC0R.NE.1)WRITE(6,AAA) 
IF<LCOR.EQ.l)WRITE(6,A33) 
WPITE(6,222) 

CALL  DRAW( 50 t X ,PR1 , 1 ,0 , LAB, IT ,0  ,0 ,0 »0 t0 , 2, 6, 12,0 t L ) 
CALL  DRAW(5J,X,PD1,2,0,LAB,IT,0,0,G,0,0,2,6,12,0,L) 
CALL  DRAW(50,X,VR1,2,0,LAB,IT,0,0,0,0,0,2,6,12,0,L) 
CALL  DRAW( 50 , X , VD1 ,3,0 ,LAB, I T ,0 ,0 ,0  ,0 ,0 ,2 ,6, 12,0,L I 
CALL  DRAW! 50, X, PR 2, 1,0, LAB, I P, 0,0, C ,0 ,0,2, 6, 12,0, L) 
CALL  DRAW(50,X,PD2,2,C,LAB, IP,0 ,0 ,0 ,0 ,0,2, 6,  12,0, L  ) 
CALL  DRAW<  5  J ,X , VD2 ,2,0 , LAB, I P,0,0 ,0 ,0 ,0,2 ,6, 12 ,0 ,L) 
CALL  DRAW(50,X,VR2,3,O,LAB, IP, 0,0,0 ,0,0, 2, 6, 12,0,L) 
CALL  DRAW(50,X,ERl,l,0,LAB,IE,O,0,0,0,O,2,6,12,O,L) 
CALL  DR  AW(  5  J,  X,  EDI,  3,0,  LAB,  IE, 0,0, CO, 0,2,  6, 12,0,  LI 
CALL  DRAW(50,X,ER2,1,0,LAB, IF ,0 ,0,0 ,0 ,0, 2, 6, 12,0, L ) 
CALL  DRAW(  5  J  ,  X  ,  ED2  ,3 ,0  ,  LAB  ,  IF  ,C  ,0  ,0  ,0  ,0  ,2  ,  6,  12,0  ,  L  ) 

30     FCRMAT(2I5) 
111   FCRMAT(9X, I2.12F10.A/) 

125    FCRMAT19X, 4F10.4/) 
222   FCPMAT(lHl) 
16C   FCRMAT(//9X, 'GAINS  FCR  FULL  ANC  RECUCED  FILTERS*//) 

161  FCRMAT(//9X, 'THEORETICAL  VARIANCE  OF  ESTIMATICN  ERROR* 
1//) 

162  F0RMAT(//9X,'REAL  VARIANCE  OF  ESTIMATION  ERRCR*//) 

163  FCRMAT<//9X, "MEANS  OF  ESTIMATION  ERROR*//) 

164  FCRMAT<//9X,*PERF0RPANCE  DEGRADATION*//) 
13  8   FORMAT (///9X,* INITIAL  I ZATI0N-E0/-1*/) 

555    FCRMAT(///9X,*K*,5X,*G<1, 1) * , 4X, *G< 2.1 ) • , 4X, *G< 3, 1) • 

1,AX,*G(A,1)*,3X,*GR(1,1)* ,3X, *GRC2, 1 1  •// ) 
666    FCRMAT(///9X,*K',4X,'PK( 1,1)*,3X,*PKC2,2)*,3X, 

2'PK(3,3)',3X,'PM4,4)',3X,'PR(1,1)*,3X,*PR<2,2)*//) 
777    FCRMAT(///9X,*K' , AX , • VK ( 1 ,1 ) * ,3X , • VK(2,2 ) * . 3X, 

2'VM3,3)',3X,*VK(A,A)*,3X,*VR< 1 , 1) • ,3X, • VR ( 2 ,2) * //) 
999    F0RMAT(///9X,*K* ,5X ,* MEANE1 * , AX , *MEANE2* . AX, *MEANE3* 

1,4X, *MEANEA*,3X,«MEANE1R* , 3X,  • MEANE2R* // ) 
99  9         FCRMAT(///9X,*K*,5X,*D(1,1)*,AX,*D(2,2)  *,AX,*E<1,1)* 

1,AX,*E<2,2)*//) 
AAA         FCRMAT1//9X, 'CORRECTION    DUE    TO    DEGRADATION    APPLIED*) 
A33         FCRMAT(//9X,'N0    CORRECTION    FOR    DEGRADATION*     I 
133         FCRMAT(9X, *PHI     MATRIX*/) 
121  F0RMAT(///9X, 'GAMMA    MATRIX') 

133  FCRMAT(///9X,' INITIAL IZAT I0N-X0/-1 •/) 

134  FORMAT <///9X,» I NI TI AL I Z ATION-PO /-l • / ) 

135  FCRMAT1///9X, 'VARIANCE    OF    FORCING    NOISE=0.01*» 

136  FCRMAT1///9X, 'VARIANCE    CF    MEASUREMENT    NOISE=0.25*) 
132         F0RMAT(///9X,'Q    MATRIX'/) 

STOP 
END 


114 


SUBROUTINE  UNIT(UtN) 
CIKENSION  U(4,4) 
CC  1  1=1, N 

CC  1  J  =  1,N 

u( i,j)=r.o 

IF(J.EQ.I)U( I, J>=1.? 

CONTINUE 

RETURN 

ENC 


SUBROUTINE  GAUSS ( I  X  ,  S, AP , V) 

A=C.C 

DO  53  1=1,12 
CALL  RANDU( IX, IY,Y) 
IX  =  IY 
5?  A  =  A+Y 

V=(A-6.C)*S+AM 

RETURN 

END 


SUBROUTINE    R/>NDU  (  I  X  ,  IY  ,  YFL  ) 

IY=IX*6553<5 
IF(IY)5,6,6 

5  IY=IY+2147483647+1 

6  YFL=IY 
YFL=YFL*.4656613E-<3 
RETURN 

ENC 


SUBROUTINE  LCC (  I , J ,  IR ,N ,M ,MS ) 

IX  =  I 
JX=J 
IF(MS-l)  1C,20,30 

io  ipx=n*( jx-n  +  ix 

GO  TO  36 
20  IFUX-JX)  22,24,24 
22  IRX=IX+( JX*JX-JX)/2 

GC  TO  36 
24  IPX=JX+( IX*lX-lX)/2 

GC  TO  36 
3 J  IPX=C 

IF(IX-JX)  36,32,36 
32  IPX=IX 
36  IP=IRX 

RETURN 

END 


SLBRCUTINE  SMPY ( A,C  ,R  ,N  ,H  ,MS ) 

CINENSICN  A(l)  ,R(1) 

CALL  LCCCN,*, IT,N,M,MS) 

DC  1  1  =  1. IT 

R(I)=A(  I)*C 

RETURN 

END 
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SIBROUTINE    MINV(A. N,D.L?M) 

CINENSION    A(l) ,L(1) ,M( ll 

C  =  1.C 

NK=-N 

CC    8f     K=1,N 

t\K  =  NK+N 

L(K)=K 

KK=NK+K 

BIGA=A(KK) 

DC    2^    J=K,N 

IZ=N*< J-l) 

CC    2T    I=K,N 

IJ=IZ+I 
10    IF(    *BS(BIGM-    ABS(A(IJ)))    15,20,20 
15    BIGA=A( I Jl 

L(K)=I 

MK)=J 
2.1    CCMINUE 

J=L(K) 

IF(J-K)  35,35,25 
25    KI=K-N 

CC    3t     1=1, N 

KI=KI+N 

FCLD=-A(KI  ) 

JI=KI-K+J 

MKI  )  =  A(JI  ) 
21    A( JI )    =HCLD 

35  I=M(K) 

IF(I-K)    45,45,38 

36  JP=N*(I-1> 
CC    AC    J=1,N 
JK=NK*J 
JI=JP+J 
HCLD=-A(JK  ) 
A< JK)=A<JI  ) 

40    A( JI )    =HOLD 

45  IF(BIGA)    48,46,48 

46  C=0.0 
RFTURN 

48    CC    55    1=1, N 

IFU-K)  50,55,50 
50     IK=NK+I 

A(  IK)  =  A(IKI/(-BIGA) 
55    CCMINUE 

CC    65    1=1, N 

IK=NK+I 

HCLD=A( IK) 

IJ=I-N 

CO    65    J=1,N 

IJ=IJ+N 

IFU-K)  60,65,60 
60  IF(J-K)  62,65,62 
62    KJ=IJ-I+K 

A(  IJ)=HCLD*A<KJ)+A(IJ) 
65    CCMINUE 

KJ=K-N 

CC    75    J=1,N 

KJ=KJ+N 

IF(J-K)  70,75,70 
7C  A(KJ)=A(KJ)/BIGA 
75    CCMINUE 

C=C*BIGA 

A(KK)=1.C/BIGA 
80    CCMINUE 

K=N 
UO    K=<K-1) 

IF(K)    150,150,105 
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1C5    I=L(K) 

IF(I-K)     12r,12C,lC8 
1C8    JC=N*(K-1) 

JR=N*< 1-1) 

DO    IK     J  =  1,N 

JK=JQ+J 

HCLD=A( JK) 

JI=JP+J 

MJK)  =  -A(J  I  ) 
110    MJI)    =HCLC 
120    J=M(K) 

IF(J-K)     It C, 100 f 12 5 
125    KI=K-N 

CO    120    1=1, N 

KI=KI+N 

HCLD=A(KI) 

JI=KI-K+J 

A(KI)=-A(JI) 
13C    MJI)    =HCLC 

GC  TO  IOC 
150  RETURN 

END 


SUBROUTINE  GMTRM  A  t  R,  N,  M  ) 
CINENSICN  Ml)  ,RU) 

IR=0 

CC    It     1=1, N 
IJ=I-N 
CO    1C     J=1,M 
IJ=IJ+N 
IR=IR+1 
1C    R(  IR)=A( IJ) 
RETURN 
ENC 


SUBROUTINE  GMPRC ( A, B, R, N, M, L ) 
DIMENSION  4(11  ,E(1)  ,R<1> 

IR=0 
IK=-M 

CO  1C  K=1,L 
IK=IK+M 
DC  10  J=1,N 
IR=IR+1 
JI=J-N 
IB  =  IK 
R(IR)=0 
CC  10  1=1, V 
JI=JI+N 
IB=IB+1 
10  R(IR)=R<IP)+A(JI)*fi<IB) 
RETURN 
ENC 


StBRCUTINE  GMSUE(A,B,R,N,M) 
CIMENSION  A(l)  ,B<1)  ,R<1) 
NF=N*M 

DC  10  I=1,NM 
10  R(  I)=A(  I)-B( I) 
RETURN 
END 
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10 


SteRCUTINE  GMADC!A,B,R,N,M) 
CIMENSION  A(l)  ,6<U  ,R(1) 

CC  1^  1=1, NM 
R(I  )=A(  I  )+B< I) 

RETURN 
END 
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APPENDIX  C 


21 

C 
C 

c 


IMPLEMENTATION  OF  A  P.ECOCEl)  SUP-OPTIMAL  FILTER  USING 
VATHIX  PSFUOO  INVERSION.  FOUR  ONE  DIMENSIONAL  FILT- 
ERS ARE  RUN  IN  PARALLEL  IN  CROC  R  TO  ESTIMATE  ONLY  TWO 
STATES  OF  A  FOURTH  ORDER  SYSTEM.  A  MONTE  CARLO  SIMULA 
TICN   IS  PERFORMED  IN  ORnER  TO  OBTAIN  APPROX I  MAT E n 
VALUES  FOP  THE  MEANS  AND   EXPERIMENTAL  COVA-UANCL'S 
OF  ESTIMATION  ERROR  FOR  THE  COMPLETE  AND  REDUCED  PRO- 
CESSES. 


CI^E 
II. HT 

2)  tVW 
3S1  (4 
45(4, 
52(4, 
6  (  1 1 4 
7,4), 
8,SE3 
<;SV33 
1LI  (2 
2Y4(5 
3,PD2 
4'  ),E 

REAL 
1  LAR 

REAL 
1  LAP 

REAL 
1  LAR 

REAL 
1  LAR 

REAL 

DATA 
ICC  5 

DATA 
2A(2, 

3)  ,A( 
4  .,- 

DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
1  .(  , 
DATA 
CATA 
DATA 
DATA 
DATA 
DATA 
CATA 
DATA 
WRIT 
CC  3 
SV33 
SV44 
SE3( 
SF4( 
SV11 
SV22 
SVR1 
SVR2 
SEK 
SE2( 
SE1P 
SE2R 
CENT 


NSICN 
(4,1). 

R(l  ,1) 
,4)  ,S2 
1),AI1 
1) ,AY3 
),A3( 1 
GH<  1,4 
(5  )  ,S 
(5  ),S 
),MM(2 
,2)  ,Y 
(5  ),V 
Dl(5i  ) 
*8  IT( 
A  BOX 
*8  IP( 
A  BOX 
*8  IE( 
A  BOX 
*8  IF( 
A  OCX 
LAB/4 
NT.NS 
/ 

A(  1,1 

4  )  ,  A  (  3 

4 , 4  )  A 

118.,- 

B  (  1  , 1 

H(  1,1 

HT(  1, 

VW(1, 

UK,SV 

VI  ,V2 


PHI  ( 

Rlli 
,  AX1 
(4,4 

(4,1 

,  H )  , 
)  ,DX 
R*(5 
V44( 
),LL 
5  (  bj 
Dl(  5 
,ED2 
12)/ 
11 


)/ 
)/ 


4,h)  , GAM (4,1)  ,A(4,4),B(4fl),Q(4,4)fH(l,4 
l),G(4,l),PKl(4,4),PK(4,4),Vw(l,i),U(4,4 
(4,4),AX2(4,4),AX3(4f4),AT(4,4),S(4,^), 
), 8X1(4,1;, 8X2(4, 1), 8X3(4,1), 8X4 (4,  1),BX 
),AI2(4,1),AI3(4,1),AI4(4,1),AY1(4,1),AY 
) ,AY4(4, 1»,AY5<4, 1),GAMU(4, l),Al(l,4),A2 
A4(l ,4) ,811 (1 .4) ,B12 (1,4) ,B13(1,4) ,B14( 1 
2(l,4),CX2< 1, 1) ,CX3( 1,1) ,SE1(5  I  ,SF2( 5  ) 

I.SE1R15  ),SF2R(5  )  ,  SVl  1  (  5.  ) ,  S  V22  (  t>  ), 
5  ).SVRll(5  I  ,SVR22(5  <)  .XAXIS(5  ),XE(4), 
L( 1  ) , MMM(  I)  ,Y1( 3  ,4) , Y2( 5  ,4) , Y3( 5  ,3) , 
,2),Y6(5  ,1),Y8(5  ,6)  ,Y9(5  ,3),Y11(5  ,1) 
.  )  ,VD2(  5  )  .SOU  5  > ,SD2<5>) ,VR1( 5  ),VR2<5 
(5  ),ER1(5  ),ER2(5  ),PD1(5'  )  ,X(  5  ) 
•VARIANCES  OF  ESTIMATION  EPRCR  FOR  XI  II 

•/ 
ION  ERR 

•/ 
RROR  FO 

•/ 
•MEAN  OF  ESTIMATION  FRROR  FO 

•/ 


•VARIANCES    OF    ^ST 
'MEAN    OF    ESTIMATI 


IMAT 

CN    El 


1 

11 

12 

11 

12)/ 

11 

H  / 

,NE, M,MR,N,NR,NZ, I  X  ,  T /4    ,5,2'     ,1,1  ,4,1  ,1,3 


OR    FOR    X(  2) 
R    X  (  1 ) 
R    X(2) 


)  ,A( 
,1), 

.,1 
19./ 
),B( 
)  ,H( 

1),H 
1  )  ,R 
1,SW 
,V3, 


1,2)  ,A( 1,3)  ,A(  1,4 
A  {  3 ,  2  )  ,  A  (  3 ,  3  I  ,  A  (  3 


•  , 


AMI,  1 


,1) 
,1) 
,1) 


,    • 

1),A 
A2(l, 1),A 
A3ll.ll  ,A 
A4(l,  1),A 
AIK1  ,1), 
A  I  2(  1 
AI3(1 
AI4(1 

E( 6,22  2) 

1    K=1,NS 

(  K  )  =    . 

(K  )=    . 

K)=i  .f 

K)=:.< 

( K  )  =    • 

(K)=-  .' 

1 ( K  )  =0  . 

2(K)=     . 

K)=.  .< 

K)=    mi, 

(K  )=    .. 

<  K  )  =*  . 

IMUE 


2,1), 
1,2)  , 
T(2,  1 
(1,1  ) 

1.AW1 
V4,AM 

/ 
K  1. 
2(1, 
3(1, 
4(1, 
AIK2 
AI2(2 
AI3(  2 
AI4(2 


B(3,  1)  ,B 
H(l,3) ,H 
),HT(3,1 
iVWRll.l 
/  .5,  .5 
1,AM2,AM 


1. 

(4,1 
(1,4 
),HT 
),KA 

,     .1 

3,  AM 


)  ,A(2,1)  ,A(2.2)  ,A(2,3)  , 
,4) ,A(4,1) ,A(4,2) ,A(4,3 

•  »-•»'   •  » !   »,l»,-4       .,-3? 


)  /     •  »     •  ,  :   .  ,  4       .  / 
)/l  .  »     •  ,     .  ,     ./ 

(  4  ,  1  )  /  1 .  ,     .  ,     .  ,     .  / 
DA/.  >1,  '.25,     .     1,1/ 
,    .     / 

4/    . 7,7,1.41, 1.41,2. 


)  ,A1(  1,3 
),A2(  1,3 
)  ,  A  3  ( 1  ,  3 


)  ,A1 
),  A2 
)  ,A3 
)  ,A4( 1,3) ,A4 
3.1) 
3,1) 
3,1) 
3,1) 


,1 ) ,AI1( 
,1)  ,AI?( 
,  1  )  ,  A  I  3  ( 
,3  )  ,AI4( 


(1,4)/!.,  . , 
(1,4) A .,1., 
(1.4)/  .,  .,1 
(1,4)/  !..?  .,  , 
,AI1(4,1)/1., 
,AI2(4,1)/,  .  , 
♦  A  I  3  (  4  ,  1  )  /  '  .  , 
,  A  I  4  (  4  ,  1  )  /     .  , 


./ 
./ 
./ 
./ 
•  , 

1.  , 


./ 
./ 
./ 

1./ 


MONTE    CARLO    SIMULATION 


119 


66 


C 
C 

c 


81 

6 


123 


54 


CC  1 
CALL 
CALL 
CALL 
CALL 
X1K  = 
X2K  = 
X3K= 
X4K  = 


L=1,NE 

GAUSS' 

GAUSS 

GAUSS 

GAUSS 

XI 

X2 

X? 

X4 


GENERATING  PHI  AND  GAM^A 


IF(L 

CALL 
UR=1 
CALL 
CALL 
CALL 
CALL 
CALL 
CC  6 
XL  =  I 
RP=1 
RP1  = 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CCNT 
WRIT 
DC  A 
WRIT 
WRIT 
WRIT 


.GT.l  ) 
UNIT( 

.( 
UNIT( 
UNIT( 
SMPY( 
SMPY( 
SMPY( 

6  1  =  1, 


GO  TO  6 

U,N) 

PHI,N) 
AX1,N) 

B,T,GAM,N,M, ;) 
B,T,BX2,N,M,f*  ) 
A,T,  ATtNtNfM 
NT 


,AT,A,N,N,N) 


,AX1,PHI,N,N) 
BX2,B,N,N,M) 


.(  /XL 
l.J/MX 

GMPRD 

SMPY(  A,RP,AX1,N,N,0) 

GMADD 

GMPRD 

SMPY( B,RP1,BX2,N,M,  ) 

GMADD ( GAM , EX 2 , GAM, N, M ) 
IMUE 
E(6,13 
C  1=1, 
E(6,12 
E(*,13 
E(6,12 


L  +  l. 
(  AX1 


(PHI 
(AT, 


) 


N 

5)  (PHK  I 

1) 

5)GAM 


J) ,J  =  1  ,N) 


»N 
tN 
.  3 


GENERATING  Q 

CALL  GMPRD(GAM,VW,BX1 ,N,M,M) 

CALL  GMTRA(GAM,DX2,N,M) 

CALL  GMPRD(BX1 ,DX2,GtN»^tNI 

WPITE(6,132) 

CC  81  1=1. N 

WRITE(6,125) (Q( I, J) ,J=l,N) 

WPITE(6,135) 

WRITE(6,136) 

CCNTINUE 

DC  123  1=1 

DT  123  J=l 

PKK  I,J)=l 

CCNTINUE 

PK1(1,1  )=C  .5 

PK  1(2, 2) =2. ) 

PK1(3,3)=2.  > 

PK1(A,4)=4. 

DC  54  1=1, N 

CC  54  J=1,N 

S(I, J)=PK1(I,J) 

XE(1  )=0.( 

XF(2>=  .< 

XE(  3)  =  !  .« 

XE(4)=c . 

INITIALIZATION  VALUES  FDR  REDUCED  FILTERS. 

C£LL  GMPRD(A1,PK1,CX2,NR,N,N) 
CALL  GMPRD (0X2 f All ,PK11 , NR.N.NR I 

CALL  GMPRD(A2,PK1,DX2,NR,N,NI 
CALL  GMPRD (DX2,AI2,PKl2,NR,N,NR) 
CALL  GMPRDU1,  XE , YE  1 , NR ,N , NR ) 
CALL  GMPRD(A2,XE,YE2,NR,N,NRI 


120 


CALL  GMPRD( A4,PK1,CX2,NR,N,N) 

CALL  GMPRD(DX2,AI4,PK14,NR,N,NR) 

CALL  GMPPD(A3,PK1,CX2,NR,N,N) 

CALL  GMPRC<DX2,AI3,PK13,NR,N,NR) 

CALL  GMPRD( A3, XE,YE3,NR,N,NR) 

CALL  GMPPD(A4,XE,YE4,NR,N,NR) 

IF(L.GT.1)G0  TO  65 

fcRITE<6,134) 

DO    142     1=1 fN 
142       WPITE(6,125HPK1(I  ,J),J=1,N) 

HPITE(6,133) 

WRITE<6,125)  XE 
65  CCNTINUE 

DC    3    K=1,NS 

CALL    GAUSS  IIXvSVltO.<  ,V1) 
CALL   GAUSSUXtSfcltAWltWl) 

IF(K.E0.1)GO    TO    4 

X1K=PHI (1,  1)*X1+PHI<1,2)*X2*PHH1,3)*X3  +  PHI(  1,4)*  X<* 
1+GAM(1,1)*(UK+W1) 

X2K=PHI<2,  lJ*Xl*PHn2t2l*X2+PHIC2t3l*X3+PHI(2t4>*X4 
2+GAM(2,l)*(UK+Wl) 

X3K=PHI<3, 1)*XH-PHI<3, 2)*X2*PHI (3 »3)*X3+PHI (3 t 4)*X4 
3+GAM(3, 1)*<UK+W1) 

X4K=PHI  (4,  1)*XH-PHI  (4,2)*X2+PHI  (4,  3  )*X3+PH  I  (  4  ,  A  )  *X4 
4  +  GAMU,  l)*(UK-i-Wl) 

MEASUREMENT 

Z=H(1,1)*X1K+V1 

GAIN  AND  COVARIANCES  OF  EST,  ERROR  FOR  COMPLETE  SYSTEM 


CALL  GMPRD(FK1,HT,BX1,N,N,NZ) 
CALL  GMPRD(H,BX1,CX2,NZ,N,NZ) 
CX2(1,1)=1.0/(CX2( 1, 1)+R(1.  1)1 

CALL  GMPRD(BX1 ,CX2,G,N,NZ,NZ) 
CALL  GMPRD(G,H,AX1,N,NZ,N) 
CALL  GMSUB(U,AX1,AX1,N,N) 
CALL  GMPRD(AX1,PK1,PK,N,N,N) 
CALL  GMTRA(PHI ,AX3,N,N) 
CALL  GMPRDtPK, AX3, AX1,N,N,NI 
CALL  GMPRD(PHI , AX1 , AX2 , N ,N,N) 
CALL  GMADD(AX2,Q,PK1,N,N) 

USING  KALMAN  FILTER  EQUATION  WITH  FULL  PLANT 


IF(K.GT 

CALL  GM 

GC  TO  8 

CALL  GM 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

bl=XlK- 

E2=X2K- 

E3=X3K- 

EA=X4K- 


SM 
GM 
GM 
CM 
GM 
SM 
GM 


.l)GO  TO  7 
PRD(U,XE,BX1 ,N,N,M) 

PRD(PHI  ,XE,BX1,N,N,1) 

PY(GAM,UK,GAMU,N,M,0) 

ADD(BX1,GAMU,BX1,N,M) 

PRDCGtHf AXltNtMt M) 

SUB(U,AX1,AX2,N,N) 

PRD(AX2,BX1,BX2,N,N,MI 

PY(G,Z,BX1, N,M,0> 

ADD(8X1,BX2,XE,N,1) 

XE(  1) 

XE(  2) 

XE(3) 

XE<4) 


CCMPUTING  PHI, GAMMA  AND  H  VALUES  FOR  REDUCED  FILTERS, 

IF(L.GT.1.0R.K.GT.l)GO  TO  61 
CALL  GMPRD<A1,PHI,B11,NR,N,N) 
CALL  GMPRD(BU,AI1,PHI1,NR,N,NR) 
CALL  GMPRDi Al , GAM, GAM1 ,NR ,N , NR  ) 
CALL  GMPRD(H,AU,H1,NR,N,NR» 
FT1=H1 
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c 
c 

c 


61 


73 


CALL  GMPRD 
CALL  GMPRD 
CALL  GMPRD 
CALL  GMPRD 
HT?=H2 
CALL  GMPRD 
CALL  G^PRD 
CALL  GMPRD 
CALL  GMPRD 
HT3=H3 
CALL  GMPRD 
CALL  GMPRD 
CALL  GMPRD 
CALL  GMPRD 
HT4=H4 

C  REDUCED 


( A2t PHIf B12, NR,N,N) 
(dl2,AI2,PHI2,NR,N,NR) 
( A2,GAM,GAM2,NR,N,NR> 
(H,AI2,H2,NR,N,NRi 


( A3,GAM,GAM3, NR,N,NR) 
><H,AI3,H3,NR,N,NR) 

(A4,PhI ,B14,NP,N,N) 
>(B14,AI4,PHI4,NR,N,NR) 
>(  A4,CAM,  GAMA,NR,N,NR) 

(H,AI4,H4,NR,N,NR) 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
C(  NT 


GMPRD 
GMTRA 
GMPRD 
GMPRD 
GMTRA 
GMPRD 
GMPRD 
GMTRA 
GMPRD 
GMPRD 
GMTRA 
GMPRD 
INUE 


(GAM1 
(GAM1 
(CXI  , 
(GAM2 
(GAM  2 
(CXI, 
(GAM3 
(GAM3 
(CXI  , 
(GAMA 
(GAMA 
(CXI  , 


tCX 
CX2 
,VW 
,CX 
CX2 

,vw 

fCX 

CX2 

t  vw 
,cx 

CX2 


R,CX 
2,NR 
tQl, 
R,CX 
2,NR 
tQ2, 
R,CX 
2,NR 
tC3, 
R,CX 
2,NR 
,G4, 


1,NR,MR,MR) 
,MR) 

NR,MR,NR) 
ltNRt MR, MP) 
,MR) 

NR,MR,NR) 
1  ,NR,MR ,MR) 
,MR  ) 

NR,MR,NR) 
1  ,NK,MR,MR) 
,  MR  ) 
NR,MR,NR) 


SLB-OPTIMAL    FILTERS 


CALL  GMPRD 
CALL  GMPRD 
CX2( 1,1)=1 
CALL  GMPRD 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GMPRD 
GMSUB 
GMPRD 
GMTRA 
GMPRD 
GMPRD 
GMADD 


ir(K.GT.l) 

YP1=YE1 

YP2=YE2 

YF3=YE3 

YPA=YE4 

GU1=GAM1*U 

G12=GAM2*U 

GL3=GAM3*U 

GU4=GAM4*U 

IF(K.EQ.l) 

CONTINUE 

CALL    GMPRD 


(PK1 
(HI, 
.  /( 
(CXI 
(Gl, 
(UR, 
(CXI 
(PHI 
(PR1 
(PHI 
(CX2 
GO    T 


K 
K 
K 
K 

GC 


1,HT1,CX1,NR,NR,NZ) 
CX1,CX2,NZ,NR,NZ) 
CX2(  1,1)+R(1,D) 
,CX2,G1,NR,NZ,NZ  ) 
H1,CX1,NR,NZ,NR) 
CX1,CX1,NR,NR) 
,PKll ,PR1,NR,NR,NR) 
1,CX2,NR,NR) 
,CX2,CX1,NR,NR,NR) 
1,CX1,CX2,NR,NR,NR) 
,Q1,PK11,NR,NR) 
0    73 


TO    37 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GMPRD 
GMPRD 
GMPRD 
GMADD 
GMADD 
GMADD 
GMPRD 
GMPRD 
GMPRD 
GMPRD 
GMPRD 


YP1=YP1+GU 
YP2=YP2+GU 
YF3=YP3+GU 
YP4=YP4+GU 


(All 

(AI2 

(AI3 

(AI4 

(  AY1 

(AY3 

(AY3 

(PHI 

(Al, 

(A2 

(A3 

(A4 

1 

2 

3 

A 


tYEl, 
,YE2, 
tYE3, 
»YEA, 
,AY2, 
t  AY4, 
tAYl, 
tAYl, 
AY5,Y 
AY5,Y 
AY5,Y 
AY5,Y 


AY1,N 
AY2,N 
AY3,N 
AYA,N 
AY1,N 
AY3,N 
AY1.N 
AY5,N 
P1,NR 
P2,NR 
P3,NR 
P4,NR 


,NR,NR) 

,NR,NR) 

,NR,NR) 

,NR,NR) 

,NR) 

,NR) 

,NR) 

,N,NR) 

,N,NR) 

,N,NR) 

,N,NR) 

,N,NR) 
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37 


ILTE-K 


CALL    GMPRD (G1,H, GH,NR,NR,N) 
C£LL    GMPRD (GH, AY2, CA,NR,N,NR  ) 
GZ=G1*Z 
YE1=GZ+YP1-GA 


FILTER     2 


CALL 
CALL 
CX2( 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GMPPD 

GMPRD 

It  1 )  =  1 

GMPFD 
GMPRD 
GMSUB 
GMPPD 
GMTRA 
GMPRD 
GMPRL 
GMADD 


(PK1 
(Hit 
.  /( 
(CXI 
(G2, 
(URf 
(CXI 
(PHI 
(PR2 
(PHI 
(CX2 


2,HT 
CXI, 
CX2( 
,CX2 
H2,C 
CXI, 
,PK1 
2,CX 
,CX2 
2,CX 
,Q2, 


(2,NP ,NP  ) 
:,CX1 ,NR,NR,NR  > 
(1,CX2,NR,NR,NR) 
PK12,NP,NR) 


CALL    GMPRD ( G 2 , H , GH , NR , N R , N ) 
CALL    G!*PRD(GH,  AY2 , G A , NR ,N ,NR ) 
GZ=G2*Z 
YE2=GZ+YP2-GA 


FILTER    3 


CALL 
CALL 
CX2( 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GMPRD 
GMPRD 
1,1)=1 
GMPRD 
GMPRD 
GMSUB 
GMPRD 
GMTRA 
GMPPD 
GMPRD 
GMADD 


(PK1 
(Hit 
.  >/( 
(CXI 
(G3, 
(URf 
(CXI 
(PHI 
(PR3 
(PHI 
(CX2 


3,HT3,CX1,NR,NR,N 
CXI ,CX?,NZ,NR,NZ) 
CX2(  1,1)+R(1,1) ) 
-r*?,G3,NR,NZ,NZ ) 
X1,NR,NZ,NR) 

r v  1       MD      km  i 


c 

,CX2  . 
H3,CX 


CALL    GMPRD(G3,H,GH,NR,NR,N) 
CALL    GMPRD(GH,AY2,GA,NR,N,NR) 
GZ=G3*Z 
YE3=GZ+YP3-GA 

FILTER    4 


CALL 
CALL 
CX2( 
CALL 
CALL 
CALL 
CALL 
CX2( 
CALL 
CALL 
CALL 


GMPPD 
GMPRD 

1,1)=1 
GMPRD 
GMPRD 
GMSUB 
GMPRD 

1,1 )=P 
GMPRD 
GMPPD 
GMADD 


HI4 

(PR4,CX2,CX1.NR,NR,NR) 
(PHI4,CX1,CX2,NR,NR,NR) 
(CX2fC4,PKI4,NR,NR) 


CALL  GMPRD(G4,H,GH,NR,NR,N) 
CALL  GMPRD (GH, AY2 , GA , NR , N , NR ) 
GZ=G4*Z 
YEA=GZ+YP4-GA 
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c 

c 
c 


la 


26 


SUB0P1 IMAL  ESTIMATES 

CALL  GMPRD(AI1,YE1,BX1,N,NR,NR) 

CALL  GMPRD( AI2,YE2,BX2,N,NR,NR) 

CALL  GMPRD< A  13 , YE3 , BX4 , N ,NR, NR  ) 

CALL  GMPRlM  AI4,YE4,BX5,N,NR,NR) 

CALL  GMAOD( 8X4,8X5, BX4,N,NR) 

C/LL  GMADD(3X1  ,BX2,BX3,N,NR) 

CALL  GMADD(3X4,BX3,BX3,N,NR) 

ElR=XlK-BX3(lt 1 ) 

EZR=X2K-BX3(2, 1) 

X1=X1K 

X2=X2K 

X3=X3K 

X4=X4K 

CALCULATING  THEORETICAL  VARIANCE  OF  ESTIMATION  ERROR 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
IF(K 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GMPRO 
GMPhD 
GMPRD 
GMPRD 
GMADD 
GMADD 
GMADD 
GMTRA 
GMPRD 
GMSUB 
GMTRA 
■  EO.l) 
GMPRD 
GMPRD 
GMPRD 
GMPRD 
GMPPD 
GMP3D 
GMPRD 
GMPRD 
GMADD 
GMADD 


(AH 
(AI2 
(AI3 
(AI4 
(BX3 
(8X1 
(3X3 
(3X3 
(BX3 
(U,A 
(AX2 
GO  T 
(S,A 
(PHI 
(AX2 
(A, A 
(AX2 
(A, A 
(3X3 
(BX2 
<S,S 
(S,S 


,G1,B 

,G2,B 

,G3,B 

,G4,B 

,8X4, 

,BX2, 

,8X5, 

tDX2, 

,F,AX 

XI  ,AX 

f  AXlf 

0  18 

X3,A, 

i At  St 

» St  At 

XI, St 

?Q» 

XI 

,R,BX 

,0X2, 

1,S,N 

2,S,N 


,N,NR,NR) 
,N,NR,NR) 
,N,NR,NR) 
,N,NR,NR) 


A, 
SI 


XI 
X2 
X3 
X4, 

BX5,N,NR) 

8X3,N,NR) 

BX3tN»NR) 

N,NR) 

1,N,NR,N) 

2,N,M 

N9N! 

N,N,M 

N,N,N) 

N,NtN) 

N,N,M 

NtNtN) 

,N,N,N) 

2tN,NR,NR) 

S2.N,NR,N) 

,N) 

tN) 


ENSEMBLE  MEAN  OF  PK/K,  PK/K  SUB-OPTIMAL  AND  GAINS 


IF(L.GT 

Y1(K,1)^ 

Y2(K,1)  = 

Y3(K,1T 

Y<5(K,1); 

Yfc(K,lh 

Y8(K,2h 

Yfi(K,3)= 

Y8(K,4)  = 

Y8(K,5)= 

Ye(K,6)= 

Y1(K,4|: 

Y2(K,4)  = 

Y3(K,2h 

Y9(K,2)  = 


1  )GO  TO 

PK(1,1 ) 

PK(2,2) 

PK(3,3  ) 

PK(4,4) 

G(  1,1) 

G(2,l) 

G(3,l) 

G(4,l) 

Gl 

G2 

S(  1, 

S(2, 

S(  3, 

S(4, 


26 


1) 
2) 
3) 
4) 


PEAL  COV.  OF  EST.  ERROR  FOR  FULL  AND  REDUCED  SYSTEMS 

SV1KK)=SV11(K)  +  E1**2 

SV22(K)=SV22(K )+E2**2 

SV33<K)=SV33(K)+E3**2 

SV44(K)=SV44(K)+E4**2 

SVR11(K)=SVR11(K)+E1P**2 

SVR22(K)=SVR22(K)+E2R**2 

MEAN  OF  ESTIMATION  ERROR  FOR  FULL  AND  REDUCED  PLANTS, 

SE1(K)=SE1(K»+E1 
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83 
55 

3 
1 


77 


883 


757 


535 
C 


SF 
St 
SF 
SF 
SE 
IF 
DC 
Yl 
Y2 

Y3 

yc 

YA 
Y5 
Y6 
Yl 

YA 
Y5 
Yi 

Y? 
XA 

cc 

cc 
cc 

UR 
VsH 
UP 
DC 
UP 

1Y8 
UR 
UR 
UR 
DC 
UR 
UP 
UR 
UP 
DC 
UP 

1Y2 
UP 
UP 
UP 
DC 
UP 

1»Y 
UP 
UP 
UR 

cr 

up 


2<K)=SF2 
3(K)=SF3 
A(K)=SEA 
!R(K)  =  SE 
2P(K)=SE 


(L.LT 

88  K 
(Kl,3 
(Kl,3 
(Kl,3 
(  K  1 ,  3 
<K1,1 
(Kl,l 
(Kltl 
KKlt 
(Kl,2 
(  K  1 ,  2 
(Kl,2 
(Kl,  2 
XIS(  K 
MINI) 
NTINU 
NTINU 
ITE(6 
ITE<6 
ITE(6 

3!  i 
ITF(6 
(K,6) 
ITE(6 
I T  E  (  6 
ITFC6 

77    K  = 
I T  E  (  6 
ITE(  6 
I T  E  (  6 
I T  E  (  6 

8  98  K 
ITE(6 
(K,2) 
ITE(6 
I T  E  (  6 
ITE(6 

757 


22 
16 
55 
=  1 
11 

22 
16 

66 
It 

11 
22 
16 
77 
=  1 
11 

22 
16 
99 
K  =  l 
ITE(6.11 
5(K,2) 
ITE(6,22 
ITE(6,16 
ITE(6T52 
535  K=l 
I  T  E  (  6  , 1 1 


NF 
=  1 
=  S 
=  S 
=  S 
=  S 
=  S 
=  S 

=s 
)  = 

=s 
=s 
=s 
=s 
)  = 


(K) 
(K) 
(K) 
1R( 
?«( 
)GQ 
,NS 
V}  1 
V22 
V33 
VhA 
El< 
E2( 
F3( 
SEA 
E1R 
E2R 
VR1 
VK2 
Kl 


2) 


♦  F2 

♦  E3 

♦  EA 

K)  ♦ E I R 

K>+E2R 
TO    55 

(  Kl) /NE 

(Kl )/NE 

(  KD/NE 

(  KD/NE 

Kl  »/NE 

KD/NE 

KD/NE 

(Kl )/NE 

(  KD/NE 

(KD/NE 

1  (  Kl  )  /NE-YMK1,2)**2 

2 ( Kl ) /NE-Y5(K1,2)**2 


) 
5) 

fNSf KADA 
DKtY3(K,D,Y8(K,2)fY8(K,3)fYe(K,AJTvri(Ktb), 

2) 
D 
6) 

NS,KACA  , 

DK,Y1(K,D,Y2(K,D1,Y3(K,D,YQ(K,D 

2) 

2) 

1) 

,NS,KADA 

DK,Yl(K,3),Y2(K,3)»Y3(K,i)fY9(K,3),YKK,Z), 

2) 
3) 
9) 

,NS,KADA 
DK,YA(KtDtY5(K,D,Y6(K,D,Yll(K,n,Y4(K,2> 

2) 

A) 

5) 

,NS,KADA 

DKfYl(K,4)  ,Y2(K,4)  ,Y3 ( K  ,  2 ) ,Y9 (K, 2 ) 


GRAPHICAC  OUTPUT 


DC 
PD 
PD 
VC 
VC 
SC 
SC 
VP 
VP 
ED 
ED 
ER 
EP 
X( 
CC 
CA 
C£ 
CA 


1  I  K= 
1(K)=Y1( 
2(K)=Y2( 
1(K)=Y1( 

2(K)=Y2( 

1(K) =Y1( 

2(K)=Y2( 

1 (K)=Y1( 

2(K)=Y2( 

1(K)=YA( 

2<K)=YM 

1(K)=YA( 

2(K)=Y5( 

K)=XAXIS 

NTINUF 

LL    DRAW( 

LC    DPAU< 

LL    DPAW( 


DNS 
K,D 
K,D 
K,3) 
K,3) 
K,A) 
K,A) 
K,2) 
K,2) 
K,l) 
K,  1) 
K,2) 
K,2) 
(K) 


tXtPDlili'  ,LAB, IT, 
tX  ,VD1  ,2,.  ,LAB,  IT, 
tX,SDl,2t'    ,LA3,IT, 


t  »  f»t2f6»12»  »  L  ) 
*t  tO»'»2»6»l2»  ♦  L  ) 
»»     t     ?     f 2  t  6  r 12  t     »L) 
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CALL    DRAW! 5^,X,VR1.3,'  t  LAB,  IT  ,  J  ,  U,  ,j  ,0  , 0 ,  2,  6 ,  12  ,  0 ,  L  ) 
CALL    DRAW!  53tXfPD2tlf0tLABf  IPff0t0fPff!Mi.f2-f6-ffl2f#  TL) 
CALL    DRAW!  5 ; , X , VD2 , 2, f , LAB , IP  , .  ,<  ,     ,     ,  i ,  2  ,  6, 12  ,i.  ,  L ) 
CALL    DRAW  (5  * ,X , SD2 , 2,0 , LAB,  IP,'   ,  W     ,-   ,  0  ,  2  ,  6,  12,  3,  L  ) 
CALL    DRAW<  5  >,X,VR2,3,;  ,  LAB,  IP,',',     , .    ,i   ,%2  ,  6,  12  ,  :  ,  L  ) 
CALL    DRAW( 5 J , XAX IS , EDI , 1 ,     ,LAB, IE ,^  ,\  ,0,  ; ,  ), 2, 6, 8,     ,L > 
CALL    DRAW<5;,XAXIS,ER1  ,3,     ,LAB,  IE,     ,.  ,0  ,  V  ,  J  ,  2  ,6,  8,  •:  ,L  ) 
CALL    DRAW! 5  ) , XAXI S  ,  ED2  , 1 ,     , LAB, IF,     ,:  ,0,  )  fnt2»6,8,0  ,L I 
CALL    DRAW!5.>,XAXIS,ER2,3,  SLAB,  IF, t  ,0,  > ,     , Q , 2, 6, 8  ,  )  ,L) 
3  i  F0RMAT(2I5) 

HI       FCPMAT(9X,  I2tl2Fl  U4/) 
125       FCPMAK^X,  AF10.4/) 
222       FCRMAT(lHl) 
52  5  FORMAT !///9X,'K',4X,'SR!l,l)',3X,'SR(2,2)»,3X, 

l'SR(3,3)  ',3X,'  SRU,4)'//> 
55  5       FORMAT !//9X,'K',5X,'G< 1  , 1) •  ,4X , • G! 2 , 1 ) • ♦ AX , • G ( 3, 1 ) • ,4X 

It 'G!^,1)',3X, »G1! 1,1) ',3X, «G2! 1,1)  •//) 
666       FORMAT! //QX,'K' ,4X , • PK ( 1 , 1 ) • , 3X , • PK (2 ,2 ) • , 3X , • PK ( 3, 3 ) • 

1,3X, 'PK(4, 4)'//) 
777       FCRMAT(//9X,'K',4X, »VK! 1 , 1 ) • , 3X , • VK { 2 , 2 ) • , 3X , 'VK!3t3) • 
1,3X,'VK!4,4)',3X,'VR!1,1)',3X,'VR<2,2)'//) 
999  rrRMAT(///9X,'K',5X,'MEANEl',4X,'MEANE2,,4X,'MEANE3' 

1, 4 X, •MEANER1 ,3X,'MEANE1R' ,3X , • MEANE2R' //  ) 
13     F0RMAT(9X, 'PHI  MATRIX'/) 
131   FORMAT (///9X, 'GAMMA  MATRIX'/) 

132  FORMAT! ///9X,' Q  MATRIX'/) 

134  FORMAT! /// 9  X,'  INIT I AL IZATION-P  i /-l ♦ /) 

133  FORMAT <///9X,»  I  MT  I  AL  I  ZAT  ION-XJ  /-l  •  /  ) 

135  F0RMAT(///9X,»  VARIANCE  CF  FORCING  N0ISE=C01') 

13o    FCPMAT!///9X, 'VARIANCE  OF  MEASUREMENT  NOISE=).25») 
16"   F0PMAT(//9X,JGAINS  FOR  FULL  AND  REDUCED  FILTERS'//) 

161  F0RMAK//9X, 'THEORETICAL  VARIANCE  OF  ESTIMATION  ERROR' 
It  IX,  •-  FULL  PLANT'//) 

162  F0RMAT(//9X,'REAL  VARIANCE  OF  ESTIMATION  ERROR'//) 

163  F0RMAT(//9X, 'MEANS  OF  ESTIMATION  ERROR'//) 

154   F0PMATC//9X, 'SUBOPTIMAL  VARIANCE  OF  ESTIMATION  ERROR' 
1//) 
STOP 
END 

SUBROUTINE  UNIT(U,N) 
DINENSICN  U!4,4) 
DC  1  1  =  1,  N 
DC  1  J=1,N 
U(  I,J)=0.( 
IF( J.FQ.I)U( I, J)=l. 
1    CONTINUE 
RETURN 
END 

SUBROUTINE  GAUSS ( I  X ,S, AM, V ) 

A  =  L  .' 

CC    5'     1=1,12 
CALL    PANDUUX,  IY,Y) 
IX=IY 
5U    A=A+Y 

V=!A-6.C)*S+AM 

RETURN 

END 

SUBROUTINE  RANDU ( I  X , I Y , YFL ) 
C 

IY=IX*65539 
IF! IY)5,6,6 

5  IY=IY+2147483647+1 

6  YFL=IY 
YFL=YFL*.4656613E-9 
RETURN 
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43 


6 
6  2 

65 


7 
75 


80 

1 
15 

1  3 


11 

12 

12 


15 


CC 
IF( 
IK  = 
M  I 

CCN 

cc 
ik  = 

HCL 

IJ  = 

CC 

IJ  = 

IK 

IF( 

KJ  = 

A(  I 

CCN 

KJ  = 

CC 

KJ  = 

IF< 

A(K 

CCN 

C  =  C 

MK 

CCN 

K  =  N 

K=< 

IF( 

I=L 

IF< 

JC  = 

JP  = 

CO 

JK  = 

HCL 

JI  = 

MJ 

A(  J 

J  =  M 

IF< 

KI  = 

CO 

Kl  = 

HCL 

JI  = 

MK 

A(  J 

GO 

PET 

END 


5  5  I  =  1 ,  N 

I-K)  5(  ,55,5 

NK+I 

K)=A( IK)/(-HIGA) 

TINUE 

65  1=1, N 

NK+I 

D=A( IK) 

I-N 

65  J=1,N 

I  J+N 

I-K  )  t>    ,65,6" 

J-K)  62,65,62 

I  J-I+K 

J  )=HCLD*A(KJ)*A(  IJ) 

TINUF 


K-N 
75  J 
KJ+N 
J-K) 


1,N 


J-K)  7.  ,75,7. 

J )=A(KJ  J/BIGA 

TINUE 

*8IGA 

K  )=1.  /PIGA 

TINUE 


K-l 

K) 

(K) 

I-K 

N*( 

N*< 

Hi 

JG+ 

D=A 

JR  + 

K)  = 

I  ) 

(K) 

J-K 

K-N 
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KI  + 

D=A 

KI- 

I)  = 

I  ) 

TO 

URN 


) 

15  , 15  ,1"5 

)  12'  ,12  ,1  3 

K-l) 

1-1) 

J  =  l  ,N 
J 

(  JK) 
J 

-A(J  I  ) 
=  HCLC 

)  K   ,  1(K  tl2  5 

1  =  1, N 
N 

(KI) 
K  +  J 

-A(JI) 
=  HCLC 
ICC 


SLBROUTINE  GMTRA ( A , R , N , M ) 
DINENSICN  Ml)  ,R(1  ) 

IP=? 

CC    If      1  =  1,  N 

IJ=I-N 

CC    1'    j=i,y 

IJ=IJ+N 

IR=IP+1 

R(  IP)  =  A(  IJ  ) 

RETURN 

END 


SUBROUTINE  GMPRD ( A . 8 , R , N , M, L ) 
UINFNSICN  A(l)  ,B(1  I  ,R(1  ) 

IR=( 

IK=-M 

CO  1<  K=1,L 
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EMC 


SUBROUTINE  LCC <  1 , J ,  IR , N,M,MS  ) 

IX=I 

JX  =  J 

IF(MS-l)  lv,2J,3 
1  i  IRX=N*< JX-1 )+IX 

GC  TO  36 
2^  IF(IX-JX)  22,24,24 
22  IPX=IX+< JX*JX-JX)/2 

GC  TO  3  6 
24  IRX=JX+( IX*lX-IX)/2 

GC  TC  36 
2^  IPX  =  » 

IF(IX-JX)  36,32,36 
32  IPX=IX 
36  IP=IRX 

PETURN 

ENC 


SUBROUTINE    SMPY ( A.C ,R , N,M, MS  ) 

CIPENSICN    Ad)  ,R(1> 

CALL    LOC(N,M,IT,N,M,MS) 

CC    1    1=1, IT 

R(I)=A(I)*C 

RETURN 

END 


SUBROUTINE    MINV ( A, N,D, L ,M ) 

DIPENSIGN    A(l)  ,L(1),M(1) 

D=1.C 

NK  =  -N 

CC    8C     K  =  1,N 

NK=NK+N 

L(K)=K 

MK)=K 

KK=NK+K 

BIGA=A(KK) 

CC    2(     J=K.N 

IZ=N* ( J-l) 

CO    20    I=K,N 

IJ=IZ-H 
10    IF(    ABS(BIGA)-    ABS(AdJ)))    15,20,20 
15    BIGA=A( IJ) 

L(K)=I 

M(K)=J 
20    CCNTINUE 

J=L(K) 

IFU-K)     35,35,25 
25    KI=K-N 

DC    3f     1=1, N 

KI=KI+N 

FCLD=-A(KI  ) 

JI=KI-K+J 

A(KI)=A< JI  ) 

3  :    A(JI  )    =HOLD 
35     I=K(K) 

IF(I-K)    45,45,38 
38    JP=N*(I-1) 

CC    4f     J  =  1,N 

JK=NK*J 

JI=JP+J 

HCLD=-A(JK) 

A( JK)=A< JI  ) 
40  A(JI)  =HOLD 
45     IF(BIGA)    48,46,48 

4  6    0=0. C 

PETURN 
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IK=IK+M 

DC  1'   J  =  1,N 

IP=IR+1 

JI=J-N 

IO  =  IK 

R(  IR)  = 

CC  T   1  =  1, M 

JI-JI+N 

IB=IB+1 

R<  IR)  =  R<  IR  )+A( JI )*R( IB) 

RETURN 

END 


SUBRCUTINE  GMSUE ( A . B  ,R . N 

DIMENSION  A(l)  ,B(1 I  ,R<1  ) 

NP=N*M 

cc  ir  1=1, NM 

R( I )=A( I)-B( I) 

RETURN 

END 


M) 


SUBRCUTINE  GMADD ( A , B , R , N , M ) 

DIMENSION  A(l)  ,R(1),R(  1) 

NP=N*y 

CO  1C  1  =  1, NM 

R(I)=A(I )+B( I) 

RETURN 

END 
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